LArSoft  v10_04_05
Liquid Argon Software toolkit - https://larsoft.org/
AuxDetReadoutGeom.cxx
Go to the documentation of this file.
1 
9 
10 #include "cetlib_except/exception.h"
12 
13 #include <limits>
14 
15 namespace {
16  constexpr auto invalid_index = std::numeric_limits<std::size_t>::max();
17 }
18 
19 namespace geo {
20 
22 
23  AuxDetReadoutInitializers AuxDetInitializer::init(std::vector<AuxDetGeo> const& ads) const
24  {
25  return initialize(ads);
26  }
27 
28  //----------------------------------------------------------------------------
30  : fADGeoToName{std::move(initializers.ADGeoToName)}
31  , fNameToADGeo{std::move(initializers.NameToADGeo)}
32  , fADGeoToChannelAndSV{std::move(initializers.ADGeoToChannelAndSV)}
33  {}
34 
35  //----------------------------------------------------------------------------
36  std::size_t AuxDetReadoutGeom::NearestAuxDet(Point_t const& point,
37  std::vector<AuxDetGeo> const& auxDets,
38  double tolerance,
39  bool throwIfAbsent) const
40  {
41  for (std::size_t a = 0; a < auxDets.size(); ++a) {
42  auto const localPoint = auxDets[a].toLocalCoords(point);
43  double const HalfCenterWidth = 0.5 * (auxDets[a].HalfWidth1() + auxDets[a].HalfWidth2());
44 
45  if (localPoint.Z() >= (-auxDets[a].Length() / 2 - tolerance) &&
46  localPoint.Z() <= (auxDets[a].Length() / 2 + tolerance) &&
47  localPoint.Y() >= (-auxDets[a].HalfHeight() - tolerance) &&
48  localPoint.Y() <= (auxDets[a].HalfHeight() + tolerance) &&
49  // if AuxDet a is a box, then HalfSmallWidth = HalfWidth
50  localPoint.X() >= (-HalfCenterWidth +
51  localPoint.Z() * (HalfCenterWidth - auxDets[a].HalfWidth2()) /
52  (0.5 * auxDets[a].Length()) -
53  tolerance) &&
54  localPoint.X() <= (HalfCenterWidth -
55  localPoint.Z() * (HalfCenterWidth - auxDets[a].HalfWidth2()) /
56  (0.5 * auxDets[a].Length()) +
57  tolerance))
58  return a;
59  } // for loop over AudDet a
60 
61  std::ostringstream msg;
62  msg << "Can't find AuxDet for position " << point;
63 
64  if (throwIfAbsent) { throw cet::exception("AuxDetReadoutGeom") << msg.str(); }
65 
66  mf::LogDebug("AuxDetReadoutGeom") << msg.str();
67  return invalid_index;
68  }
69 
70  // ----------------------------------------------------------------------------
72  std::vector<AuxDetGeo> const& auxDets,
73  double tolerance,
74  bool throwIfAbsent) const
75  {
76  std::size_t const auxDetIdx = NearestAuxDet(point, auxDets, tolerance, throwIfAbsent);
77  if (auxDetIdx == invalid_index) { return invalid_index; }
78 
79  AuxDetGeo const& adg = auxDets[auxDetIdx];
80  for (std::size_t a = 0; a < adg.NSensitiveVolume(); ++a) {
81  AuxDetSensitiveGeo const& adsg = adg.SensitiveVolume(a);
82  auto const localPoint = adsg.toLocalCoords(point);
83 
84  double const HalfCenterWidth = 0.5 * (adsg.HalfWidth1() + adsg.HalfWidth2());
85 
86  if (localPoint.Z() >= -(adsg.Length() / 2 + tolerance) &&
87  localPoint.Z() <= (adsg.Length() / 2 + tolerance) &&
88  localPoint.Y() >= -adsg.HalfHeight() - tolerance &&
89  localPoint.Y() <= adsg.HalfHeight() + tolerance &&
90  // if AuxDet a is a box, then HalfSmallWidth = HalfWidth
91  localPoint.X() >=
92  -HalfCenterWidth +
93  localPoint.Z() * (HalfCenterWidth - adsg.HalfWidth2()) / (0.5 * adsg.Length()) -
94  tolerance &&
95  localPoint.X() <=
96  HalfCenterWidth -
97  localPoint.Z() * (HalfCenterWidth - adsg.HalfWidth2()) / (0.5 * adsg.Length()) +
98  tolerance)
99  return a;
100  } // for loop over AuxDetSensitive a
101 
102  std::ostringstream msg;
103  msg << "Can't find AuxDetSensitive for position " << point;
104 
105  if (throwIfAbsent) { throw cet::exception("AuxDetReadoutGeom") << msg.str(); }
106 
107  mf::LogDebug("AuxDetReadoutGeom") << msg.str();
108  return invalid_index;
109  }
110 
111  //----------------------------------------------------------------------------
112  // the first member of the pair is the index in the auxDets vector for the AuxDetGeo,
113  // the second member is the index in the vector of AuxDetSensitiveGeos for that AuxDetGeo
114  std::pair<std::size_t, std::size_t> AuxDetReadoutGeom::ChannelToSensitiveAuxDet(
115  std::string const& detName,
116  std::uint32_t channel) const
117  {
118  std::size_t adGeoIdx = DetNameToAuxDet(detName);
119 
120  // look for the index of the sensitive volume for the given channel
121  auto it = fADGeoToChannelAndSV.find(adGeoIdx);
122  if (it == fADGeoToChannelAndSV.cend()) {
123  throw cet::exception("Geometry")
124  << "Given AuxDetSensitive channel, " << channel
125  << ", cannot be found in vector associated to AuxDetGeo index: " << adGeoIdx
126  << ". Vector has size " << it->second.size();
127  }
128 
129  // get the vector of channels to AuxDetSensitiveGeo index
130  if (channel < it->second.size()) return {adGeoIdx, it->second[channel].second};
131 
132  throw cet::exception("Geometry") << "Given AuxDetGeo with index " << adGeoIdx
133  << " does not correspond to any vector of sensitive volumes";
134  }
135 
136  //----------------------------------------------------------------------------
138  std::uint32_t const channel,
139  std::string const& auxDetName,
140  std::vector<geo::AuxDetGeo> const& auxDets) const
141  {
142  // Figure out which detector we are in
143  auto ad_it = fNameToADGeo.find(auxDetName);
144  if (ad_it == fNameToADGeo.cend()) {
145  throw cet::exception("CRTWireReadoutGeom") << "No AuxDetGeo with name " << auxDetName;
146  }
147  std::size_t const ad = ad_it->second;
148 
149  // Get the vector of channel and sensitive volume pairs
150  auto it = fADGeoToChannelAndSV.find(ad);
151  if (it == fADGeoToChannelAndSV.end()) {
152  throw cet::exception("CRTWireReadoutGeom") << "No entry in channel and sensitive volume"
153  << " map for AuxDet index " << ad << " bail";
154  }
155 
156  // Loop over the vector of channel and sensitive volumes to determine the sensitive
157  // volume for this channel. Then get the origin of the sensitive volume in the world
158  // coordinate system.
159  for (auto const& [ch, svindex] : it->second) {
160  if (ch == channel) { return auxDets[ad].SensitiveVolume(svindex).GetCenter(); }
161  }
162 
163  return {};
164  }
165 
166  //----------------------------------------------------------------------------
167  std::size_t AuxDetReadoutGeom::DetNameToAuxDet(std::string const& detName) const
168  {
169  // loop over the map of AuxDet names to Geo object numbers to determine which auxdet
170  // we have. If no name in the map matches the provided string, throw an exception;
171  // the list of AuxDetGeo passed as argument is ignored! Note that fADGeoToName must
172  // have been updated by a derived class.
173  auto it = std::find_if(fADGeoToName.begin(), fADGeoToName.end(), [&detName](auto const& pr) {
174  return pr.second == detName;
175  });
176  if (it == fADGeoToName.end()) {
177  throw cet::exception("Geometry") << "No AuxDetGeo matching name: " << detName;
178  }
179  return it->first;
180  }
181 }
float Length(const PFPStruct &pfp)
Definition: PFPUtils.cxx:3267
Encapsulate the geometry of the sensitive portion of an auxiliary detector .
AuxDetSensitiveGeo const & SensitiveVolume(size_t sv) const
Definition: AuxDetGeo.h:123
std::size_t NearestAuxDet(Point_t const &point, std::vector< AuxDetGeo > const &auxDets, double tolerance=0, bool throwIfAbsent=true) const
std::map< std::size_t, std::string > ADGeoToName
map the AuxDetGeo index to the name
std::map< std::size_t, std::vector< chanAndSV > > fADGeoToChannelAndSV
std::size_t DetNameToAuxDet(std::string const &detName) const
std::pair< std::size_t, std::size_t > ChannelToSensitiveAuxDet(std::string const &detName, std::uint32_t channel) const
LocalPoint_t toLocalCoords(Point_t const &world) const
Transform point from world frame to local auxiliary detector frame.
virtual ~AuxDetInitializer()
size_t NSensitiveVolume() const
Definition: AuxDetGeo.h:124
virtual AuxDetReadoutInitializers initialize(std::vector< AuxDetGeo > const &ads) const =0
Encapsulate the geometry of an auxiliary detector.
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
MaybeLogger_< ELseverityLevel::ELsev_success, false > LogDebug
std::map< std::string, std::size_t > fNameToADGeo
map the names to the AuxDetGeo index
AuxDetReadoutInitializers init(std::vector< AuxDetGeo > const &ads) const
Point_t AuxDetChannelToPosition(std::uint32_t channel, std::string const &auxDetName, std::vector< AuxDetGeo > const &auxDets) const
second_as<> second
Type of time stored in seconds, in double precision.
Definition: spacetime.h:82
ROOT libraries.
std::size_t NearestSensitiveAuxDet(Point_t const &point, std::vector< AuxDetGeo > const &auxDets, double tolerance=0, bool throwIfAbsent=true) const
std::map< std::size_t, std::string > fADGeoToName
map the AuxDetGeo index to the name
Interface to auxiliary-detector geometry for wire readouts. .
AuxDetReadoutGeom(AuxDetReadoutInitializers initializers)
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33