LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
DriftPartitions.cxx File Reference

Classes describing partition of cryostat volume. More...

#include "larcorealg/Geometry/DriftPartitions.h"
#include "larcorealg/CoreUtils/DumpUtils.h"
#include "larcorealg/CoreUtils/RealComparisons.h"
#include "larcorealg/Geometry/CryostatGeo.h"
#include "larcorealg/Geometry/Decomposer.h"
#include "larcorealg/Geometry/PlaneGeo.h"
#include "larcorealg/Geometry/TPCGeo.h"
#include "cetlib_except/exception.h"
#include "messagefacility/MessageLogger/MessageLogger.h"
#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstdlib>
#include <iterator>
#include <memory>
#include <type_traits>
#include <utility>
#include <vector>

Go to the source code of this file.

Classes

struct  TPCgroup_t
 
struct  TPCwithArea_t
 
struct  SorterByKey< Key, ExtractKey, Comparer >
 

Typedefs

using TPCandPos_t = std::pair< double, geo::TPCGeo const * >
 
template<geo::part::AreaOwner::AreaRangeMember_t Range>
using SortTPCareaByAreaRangeLower = SorterByKey< double, geo::part::details::RangeLowerBoundExtractor< Range >>
 
using SortTPCwithAreaByWidth = SortTPCareaByAreaRangeLower<&geo::part::AreaOwner::Area_t::width >
 
using SortTPCwithAreaByDepth = SortTPCareaByAreaRangeLower<&geo::part::AreaOwner::Area_t::depth >
 

Functions

std::vector< std::pair< geo::DriftPartitions::DriftDir_t, std::vector< geo::TPCGeo const * > > > groupTPCsByDriftDir (geo::CryostatGeo const &cryo)
 
std::vector< TPCandPos_tsortTPCsByDriftCoord (std::vector< geo::TPCGeo const * > const &TPCs, geo::DriftPartitions::Decomposer_t const &decomp)
 
std::vector< TPCgroup_tgroupByDriftCoordinate (std::vector< TPCandPos_t > const &TPCs)
 
unsigned int checkTPCcoords (std::vector< geo::TPCGeo const * > const &TPCs)
 
template<typename Range >
geo::DriftPartitions::DriftDir_t detectGlobalDriftDir (Range &&directions)
 
geo::part::AreaOwner::Area_t TPCarea (geo::TPCGeo const &TPC, geo::DriftPartitions::Decomposer_t const &decomposer)
 
std::vector< TPCwithArea_taddAreaToTPCs (std::vector< geo::TPCGeo const * > const &TPCs, geo::DriftPartitions::Decomposer_t const &decomposer)
 
template<typename BeginIter , typename EndIter >
geo::part::AreaOwner::Area_t computeTotalArea (BeginIter TPCbegin, EndIter TPCend)
 
template<geo::part::AreaOwner::AreaRangeMember_t sortingRange, typename BeginIter , typename EndIter >
std::vector< std::vector< TPCwithArea_t const * >::const_iteratorgroupTPCsByRangeCoord (BeginIter beginTPCwithArea, EndIter endTPCwithArea)
 
template<geo::part::AreaOwner::AreaRangeMember_t sortingRange, typename BeginIter , typename EndIter >
std::pair< std::vector< TPCwithArea_t const * >, std::vector< std::vector< TPCwithArea_t const * >::const_iterator > > sortAndGroupTPCsByRangeCoord (BeginIter beginTPCwithArea, EndIter endTPCwithArea)
 
template<typename BeginIter , typename EndIter , typename TPCendIter , typename SubpartMaker >
geo::part::Partition< geo::TPCGeo const >::Subpartitions_t createSubpartitions (BeginIter itTPCbegin, EndIter itTPCend, TPCendIter TPCend, SubpartMaker subpartMaker)
 
auto makeTPCPartitionElement (TPCwithArea_t const &TPCinfo)
 
template<typename BeginIter , typename EndIter >
std::unique_ptr< geo::part::Partition< geo::TPCGeo const > > makeWidthPartition (BeginIter beginTPCwithArea, EndIter endTPCwithArea)
 
template<typename BeginIter , typename EndIter >
std::unique_ptr< geo::part::Partition< geo::TPCGeo const > > makeDepthPartition (BeginIter beginTPCwithArea, EndIter endTPCwithArea)
 
template<typename BeginIter , typename EndIter >
std::unique_ptr< geo::part::Partition< geo::TPCGeo const > > makeGridPartition (BeginIter beginTPCwithArea, EndIter endTPCwithArea)
 
template<typename BeginIter , typename EndIter >
std::unique_ptr< geo::part::Partition< geo::TPCGeo const > > makePartition (BeginIter beginTPCwithArea, EndIter endTPCwithArea)
 
template<typename TPCPartitionResultType , geo::part::AreaOwner::AreaRangeMember_t Range, typename BeginIter , typename EndIter , typename SubpartMaker >
std::unique_ptr< geo::part::Partition< geo::TPCGeo const > > makeSortedPartition (BeginIter beginTPCwithArea, EndIter endTPCwithArea, SubpartMaker subpartMaker)
 
template<typename BeginIter , typename EndIter >
auto makeCPointerVector (BeginIter b, EndIter e)
 
template<typename T >
auto makeCPointerVector (std::vector< T > const &v)
 
std::unique_ptr< geo::DriftPartitions::TPCPartition_tmakePartition (std::vector< TPCwithArea_t > const &TPCs)
 

Detailed Description

Classes describing partition of cryostat volume.

Author
Gianluca Petrillo (petri.nosp@m.llo@.nosp@m.fnal..nosp@m.gov)
Date
July 13, 2017
See also
DriftPartitions.h

Definition in file DriftPartitions.cxx.

Typedef Documentation

template<geo::part::AreaOwner::AreaRangeMember_t Range>
using SortTPCareaByAreaRangeLower = SorterByKey<double, geo::part::details::RangeLowerBoundExtractor<Range>>

Definition at line 457 of file DriftPartitions.cxx.

using TPCandPos_t = std::pair<double, geo::TPCGeo const*>

Definition at line 147 of file DriftPartitions.cxx.

Function Documentation

std::vector<TPCwithArea_t> addAreaToTPCs ( std::vector< geo::TPCGeo const * > const &  TPCs,
geo::DriftPartitions::Decomposer_t const &  decomposer 
)

Definition at line 397 of file DriftPartitions.cxx.

References TPCarea().

Referenced by geo::buildDriftVolumes().

399 {
400  /*
401  * Transforms a collection of TPCs into a collection of TPCs with area.
402  */
403  std::vector<TPCwithArea_t> result;
404  result.reserve(TPCs.size());
405 
406  for (auto const& TPC : TPCs)
407  result.emplace_back(TPCarea(*TPC, decomposer), TPC);
408 
409  return result;
410 } // addAreaToTPCs()
geo::part::AreaOwner::Area_t TPCarea(geo::TPCGeo const &TPC, geo::DriftPartitions::Decomposer_t const &decomposer)
unsigned int checkTPCcoords ( std::vector< geo::TPCGeo const * > const &  TPCs)

Definition at line 282 of file DriftPartitions.cxx.

References geo::vect::dot(), geo::DriftPartitions::driftCoord(), geo::TPCGeo::DriftDir(), e, geo::TPCGeo::ID(), lar::util::makeVector3DComparison(), lar::util::RealComparisons< RealType >::nonEqual(), and lar::dump::vector3D().

Referenced by geo::buildDriftVolumes().

283 {
284  /*
285  * Verify coordinate system consistency between TPCs:
286  * * need to have the same drift direction
287  * * need to have the same drift coordinate
288  *
289  * On error, it prints information on the error stream ("GeometryPartitions").
290  * It returns the number of errors found.
291  */
292 
293  auto iTPC = TPCs.cbegin(), tend = TPCs.cend();
294  if (iTPC == tend) {
295  mf::LogProblem("GeometryPartitions") << "checkTPCcoords() got an empty partition.";
296  return 0;
297  }
298 
299  geo::TPCGeo const& refTPC = **iTPC;
300  decltype(auto) refDriftDir = refTPC.DriftDir();
301 
302  auto driftCoord = [&refDriftDir](geo::TPCGeo const& TPC) {
303  return geo::vect::dot(TPC.FirstPlane().GetCenter(), refDriftDir);
304  };
305 
306  auto const refDriftPos = driftCoord(refTPC);
307 
309  auto vectorIs = lar::util::makeVector3DComparison(coordIs);
310 
311  unsigned int nErrors = 0U;
312  while (++iTPC != tend) {
313  geo::TPCGeo const& TPC = **iTPC;
314 
315  if (vectorIs.nonEqual(TPC.DriftDir(), refDriftDir)) {
316  mf::LogProblem("GeometryPartitions")
317  << "Incompatible drift directions between " << TPC.ID() << " "
318  << lar::dump::vector3D(TPC.DriftDir()) << " and " << refTPC.ID() << " "
319  << lar::dump::vector3D(refTPC.DriftDir());
320  ++nErrors;
321  }
322  auto const driftPos = driftCoord(TPC);
323  if (coordIs.nonEqual(driftPos, refDriftPos)) {
324  mf::LogProblem("GeometryPartitions")
325  << "Incompatible drift coordinate between " << TPC.ID() << " (" << driftPos << "( and "
326  << refTPC.ID() << " (" << refDriftPos << ")";
327  ++nErrors;
328  }
329  } // while
330  return nErrors;
331 } // checkTPCcoords()
geo::TPCID const & ID() const
Returns the identifier of this TPC.
Definition: TPCGeo.h:270
auto vector3D(Vector3D const &v)
Returns a manipulator which will print the specified vector.
Definition: DumpUtils.h:326
Provides simple real number checks.
Geometry information for a single TPC.
Definition: TPCGeo.h:36
MaybeLogger_< ELseverityLevel::ELsev_error, true > LogProblem
auto makeVector3DComparison(RealType threshold)
Creates a Vector3DComparison from a RealComparisons object.
constexpr auto dot(Vector const &a, OtherVector const &b)
Return cross product of two vectors.
Float_t e
Definition: plot.C:35
Vector_t DriftDir() const
Returns the direction of the drift (vector pointing toward the planes).
Definition: TPCGeo.h:124
Namespace collecting geometry-related classes utilities.
template<typename BeginIter , typename EndIter >
geo::part::AreaOwner::Area_t computeTotalArea ( BeginIter  TPCbegin,
EndIter  TPCend 
)

Definition at line 414 of file DriftPartitions.cxx.

References lar::util::simple_geo::Rectangle< Data >::extendToInclude().

Referenced by makeSortedPartition().

415 {
416  /*
417  * Computes an area covering the areas from all TPCs delimited by the
418  * iterators in argument.
419  * Each iterator points to a AreaOwner pointer.
420  */
421  auto iTPC = TPCbegin;
422  geo::part::AreaOwner::Area_t totalArea((*iTPC)->area());
423  while (++iTPC != TPCend)
424  totalArea.extendToInclude((*iTPC)->area());
425  return totalArea;
426 } // computeTotalArea()
void extendToInclude(Rectangle_t const &r)
Extends the range to include the specified point.
Definition: SimpleGeo.h:545
template<typename BeginIter , typename EndIter , typename TPCendIter , typename SubpartMaker >
geo::part::Partition<geo::TPCGeo const>::Subpartitions_t createSubpartitions ( BeginIter  itTPCbegin,
EndIter  itTPCend,
TPCendIter  TPCend,
SubpartMaker  subpartMaker 
)

Definition at line 547 of file DriftPartitions.cxx.

References makeTPCPartitionElement().

Referenced by makeSortedPartition().

552 {
553  /*
554  * Internal helper to create a sequence of partitions (geo::part::Partition)
555  * from groups of TPCs. Each TPC is specified as a pointer to TPCwithArea_t
556  * object.
557  *
558  * The groups are specified in a way that is more or less convenient after
559  * calling groupTPCsByRangeCoord(): the uber-iterators itTPCbegin and itTPCend
560  * delimit a collection of iterators pointing to the first TPC of a group
561  * (the pointed TPCs are on a different collection). This can define all
562  * groups, each group delimited by the TPC pointed by an iterator pointing at
563  * the beginning of that group and the TPC pointed by the next iterator,
564  * pointing at the beginning of the next group. The exception is the last
565  * group, for which there is no "next iterator pointing to the beginning of
566  * the next group". That information is provided separately by the third
567  * iterator argument, that is an iterator of a different type than the other
568  * two (which are "uber-iterators" whose values are iterators), and which
569  * points to the last available TPC, that is also the end iterator for the
570  * TPCs of the last group.
571  *
572  */
574 
575  // TPCgroups has an iterator to the starting TPC of each group.
576  // Iterators refer to the `TPCs` collection.
577  // The end iterator of the group is the starting iterator of the next one;
578  // the last group includes all the remaining TPCs and its end iterator is
579  // the end iterator of `TPCs`.
580  auto igbegin = itTPCbegin;
581  while (igbegin != itTPCend) {
582  auto const gbegin = *igbegin;
583  auto const gend = (++igbegin == itTPCend) ? TPCend : *igbegin;
584 
585  //
586  // create a partition from the new group
587  //
588  if (std::distance(gbegin, gend) == 1) {
589  subparts.emplace_back(makeTPCPartitionElement(**gbegin));
590  }
591  else {
592  auto subpart = subpartMaker(gbegin, gend);
593  if (!subpart) return {}; // failure!!
594 
595  subparts.emplace_back(std::move(subpart));
596  }
597  } // while
598  return subparts;
599 } // createSubpartitions()
std::vector< std::unique_ptr< Partition_t const >> Subpartitions_t
Type of list of subpartitions. It needs to preserve polymorphism.
Definition: Partitions.h:191
auto makeTPCPartitionElement(TPCwithArea_t const &TPCinfo)
template<typename Range >
geo::DriftPartitions::DriftDir_t detectGlobalDriftDir ( Range &&  directions)

Definition at line 335 of file DriftPartitions.cxx.

References util::abs(), util::cbegin(), util::cend(), dir, geo::vect::dot(), e, lar::util::RealComparisons< RealType >::equal(), and lar::dump::vector3D().

Referenced by geo::buildDriftVolumes().

336 {
337  /*
338  * Returns the first of the specified directions (possibly flipped);
339  * throws an exception if any of them is not parallel to that one.
340  */
341  using std::cbegin;
342  using std::cend;
343  auto iDir = cbegin(directions);
344  auto dend = cend(directions);
345  if (!(iDir != dend)) {
346  throw cet::exception("buildDriftVolumes") << "detectGlobalDriftDir(): no TPCs provided!\n";
347  }
348 
350  auto compatibleDir = [comp](auto const& a, auto const& b) {
351  return comp.equal(std::abs(geo::vect::dot(a, b)), +1.0);
352  };
353 
354  auto const dir = *(iDir++);
355  for (; iDir != dend; ++iDir) {
356  if (compatibleDir(dir, *iDir)) continue;
357  throw cet::exception("buildDriftVolumes")
358  << "Found drift directions not compatible: " << lar::dump::vector3D(dir) << " and "
359  << lar::dump::vector3D(*iDir) << "\n";
360  } // for
361 
362  // mildly prefer positive directions
363  return ((dir.X() <= 0.0) && (dir.Y() <= 0.0) && (dir.Z() <= 0.0)) ? -dir : dir;
364 } // detectGlobalDriftDir()
auto vector3D(Vector3D const &v)
Returns a manipulator which will print the specified vector.
Definition: DumpUtils.h:326
decltype(auto) constexpr cend(T &&obj)
ADL-aware version of std::cend.
Definition: StdUtils.h:93
Provides simple real number checks.
constexpr auto abs(T v)
Returns the absolute value of the argument.
constexpr auto dot(Vector const &a, OtherVector const &b)
Return cross product of two vectors.
TDirectory * dir
Definition: macro.C:5
decltype(auto) constexpr cbegin(T &&obj)
ADL-aware version of std::cbegin.
Definition: StdUtils.h:85
Float_t e
Definition: plot.C:35
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
std::vector<TPCgroup_t> groupByDriftCoordinate ( std::vector< TPCandPos_t > const &  TPCs)

Definition at line 241 of file DriftPartitions.cxx.

References geo::TPCGeo::Nplanes(), and geo::TPCGeo::Plane0Pitch().

Referenced by geo::buildDriftVolumes().

242 {
243  /*
244  * Produces a list of TPC groups. Within each group, the TPCs have similar
245  * drift coordinate.
246  * Similar is defined arbitrarily as ten times the plane pitch (of the first
247  * planes of the TPC).
248  */
249  if (TPCs.empty()) return {};
250 
251  geo::TPCGeo const& firstTPC = *(TPCs.front().second);
252  // arbitrary 5 cm if the first TPC has only one plane (pixel readout?);
253  // protect against the case where planes have the same position
254  // (e.g. dual phase)
255  double const groupThickness =
256  10.0 * std::min(((firstTPC.Nplanes() > 1) ? firstTPC.Plane0Pitch(1) : 0.5), 0.1);
257 
258  auto iFirstTPC = TPCs.cbegin(), tend = TPCs.cend();
259 
260  std::vector<TPCgroup_t> result;
261  while (iFirstTPC != tend) {
262  double const posEnd = iFirstTPC->first + groupThickness; // not beyond here
263  double sumPos = 0.0;
264  std::vector<geo::TPCGeo const*> TPCs;
265  auto iEndGroup = iFirstTPC;
266  do {
267  TPCs.push_back(iEndGroup->second);
268  sumPos += iEndGroup->first;
269  ++iEndGroup;
270  } while ((iEndGroup != tend) && (iEndGroup->first < posEnd));
271 
272  double const averagePos = sumPos / TPCs.size();
273  result.emplace_back(averagePos, std::move(TPCs));
274 
275  iFirstTPC = iEndGroup;
276  } // while (outer)
277 
278  return result;
279 } // groupByDriftCoordinate()
unsigned int Nplanes() const
Number of planes in this tpc.
Definition: TPCGeo.h:137
Geometry information for a single TPC.
Definition: TPCGeo.h:36
double Plane0Pitch(unsigned int p) const
Returns the center of the TPC volume in world coordinates [cm].
Definition: TPCGeo.cxx:314
std::vector<std::pair<geo::DriftPartitions::DriftDir_t, std::vector<geo::TPCGeo const*> > > groupTPCsByDriftDir ( geo::CryostatGeo const &  cryo)

Definition at line 168 of file DriftPartitions.cxx.

References geo::TPCGeo::DriftDir(), e, lar::util::makeVector3DComparison(), geo::CryostatGeo::NTPC(), geo::vect::rounded01(), lar::util::RealComparisons< RealType >::threshold, and geo::CryostatGeo::TPC().

Referenced by geo::buildDriftVolumes().

169 {
170  /*
171  * TPCs are grouped by their drift direction, allowing for a small rounding
172  * error.
173  * The result is a collection with one element for each group of TPCs sharing
174  * the same drift direction. Each element is a pair with that drift direction
175  * first, and a collection of all TPCs with that drift direction.
176  * Elements are in no particular order.
177  */
178  std::vector<std::pair<geo::DriftPartitions::DriftDir_t, std::vector<geo::TPCGeo const*>>> result;
179 
181  auto vectorIs = lar::util::makeVector3DComparison(coordIs);
182 
183  auto const nTPCs = cryo.NTPC();
184  for (unsigned int iTPC = 0; iTPC < nTPCs; ++iTPC) {
185  geo::TPCGeo const& TPC = cryo.TPC(iTPC);
186 
187  decltype(auto) driftDir = TPC.DriftDir();
188 
189  std::size_t iGroup = 0;
190  for (; iGroup < result.size(); ++iGroup) {
191  if (vectorIs.nonEqual(driftDir, result[iGroup].first)) continue;
192  result[iGroup].second.push_back(&TPC);
193  break;
194  } // for
195 
196  // if we did not find a group yet, make a new one
197  if (iGroup == result.size()) {
198  result.emplace_back(geo::vect::rounded01(driftDir, coordIs.threshold),
199  std::vector<geo::TPCGeo const*>{&TPC});
200  } // if
201 
202  } // for
203 
204  return result;
205 } // groupTPCsByDriftDir()
Provides simple real number checks.
Geometry information for a single TPC.
Definition: TPCGeo.h:36
for(Int_t i=0;i< nentries;i++)
Definition: comparison.C:30
STL namespace.
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:101
auto makeVector3DComparison(RealType threshold)
Creates a Vector3DComparison from a RealComparisons object.
Vector rounded01(Vector const &v, Scalar tol)
Returns a vector with all components rounded if close to 0, -1 or +1.
Float_t e
Definition: plot.C:35
template<geo::part::AreaOwner::AreaRangeMember_t sortingRange, typename BeginIter , typename EndIter >
std::vector<std::vector<TPCwithArea_t const*>::const_iterator> groupTPCsByRangeCoord ( BeginIter  beginTPCwithArea,
EndIter  endTPCwithArea 
)

Definition at line 466 of file DriftPartitions.cxx.

References lar::util::RealComparisons< RealType >::nonSmaller().

469 {
470  /*
471  * Groups each TPC with all the following ones overlapping in the selected
472  * sorting range.
473  * The range of TPCs must be already sorted by lower sorting range coordinate.
474  * The result is a list of iterators like the ones in input (in fact, the
475  * first is always beginTPCwithArea).
476  * The iterators are expected to be valid also after this function has
477  * returned (lest the result be unusable).
478  */
479 
481 
482  // tolerate 1mm overlaps; this is way forgiving, but apparently there are
483  // geometries around (DUNE 35t) with overlaps of that size.
485 
486  auto gbegin = beginTPCwithArea;
487  while (gbegin != endTPCwithArea) {
488 
489  groupStart.push_back(gbegin);
490 
491  //
492  // detect the end of this group
493  //
494  auto range = (*gbegin)->area().*sortingRange;
495  auto gend = gbegin;
496  while (++gend != endTPCwithArea) {
497  //
498  // check if the sorting range of this TPC (gend) overlaps the accumulated
499  // one; since TPCs are sorted by lower end of the range, gend has that one
500  // larger than the accumulated one, and overlap happens only if that lower
501  // bound is smaller than the upper bound of the accumulated range;
502  // we need to avoid rounding errors: close borders are decided as
503  // non-overlapping
504  //
505  auto const& TPCrange = (*gend)->area().*sortingRange;
506  if (coordIs.nonSmaller(TPCrange.lower, range.upper)) break;
507  range.extendToInclude(TPCrange);
508  } // while (inner)
509 
510  // prepare for the next group
511  gbegin = gend;
512  } // while (outer)
513 
514  return groupStart;
515 
516 } // groupTPCsByRangeCoord<>()
Provides simple real number checks.
intermediate_table::const_iterator const_iterator
template<typename BeginIter , typename EndIter >
auto makeCPointerVector ( BeginIter  b,
EndIter  e 
)

Definition at line 1001 of file DriftPartitions.cxx.

Referenced by makeCPointerVector(), and makePartition().

1002 {
1003  using value_type = typename BeginIter::value_type;
1004  std::vector<value_type const*> result;
1005  result.reserve(std::distance(b, e));
1006  std::transform(b, e, std::back_inserter(result), [](auto& obj) { return std::addressof(obj); });
1007  return result;
1008 } // makeCPointerVector()
Float_t e
Definition: plot.C:35
template<typename T >
auto makeCPointerVector ( std::vector< T > const &  v)

Definition at line 1011 of file DriftPartitions.cxx.

References makeCPointerVector().

1012 {
1013  return makeCPointerVector(v.cbegin(), v.cend());
1014 }
auto makeCPointerVector(BeginIter b, EndIter e)
template<typename BeginIter , typename EndIter >
std::unique_ptr< geo::part::Partition< geo::TPCGeo const > > makeDepthPartition ( BeginIter  beginTPCwithArea,
EndIter  endTPCwithArea 
)

Definition at line 698 of file DriftPartitions.cxx.

References lar::util::simple_geo::Rectangle< double >::depth.

Referenced by makePartition(), and makeTPCPartitionElement().

701 {
702  return makeSortedPartition<geo::part::DepthPartition<geo::TPCGeo const>,
704  beginTPCwithArea, endTPCwithArea, &makeWidthPartition<BeginIter, EndIter>);
705 } // makeDepthPartition()
Range_t depth
Range along depth direction.
Definition: SimpleGeo.h:430
template<typename BeginIter , typename EndIter >
std::unique_ptr< geo::part::Partition< geo::TPCGeo const > > makeGridPartition ( BeginIter  beginTPCwithArea,
EndIter  endTPCwithArea 
)

Definition at line 709 of file DriftPartitions.cxx.

References if(), makePartition(), lar::util::RealComparisons< RealType >::nonSmaller(), part, r, lar::util::RealComparisons< RealType >::strictlyGreater(), and lar::util::RealComparisons< RealType >::strictlySmaller().

Referenced by makePartition(), and makeTPCPartitionElement().

712 {
713  /*
714  * Requires at least 4 input TPCs (otherwise, do not use GridPartition).
715  *
716  * 1. attempt a partition on width
717  * 1. if failed, return failure
718  * 2. attempt to partition the first subpartition on depth
719  * 1. if failed, return failure
720  * 3. extend each depth partition to the other width partitions
721  * 1. if the extension of one partition line fails, discard it
722  * 2. if no partition line survives, return failure
723  * 4. run makePartition() on each of the cells
724  * 5. create and return a GridPartition object from the subpartitions so
725  * created
726  *
727  * This algorithm could use some factorization...
728  */
729  using Area_t = geo::part::AreaOwner::Area_t;
730 
731  //
732  // sort by width coordinate; work with pointers for convenience
733  //
734  auto const TPCgroupInfo =
735  sortAndGroupTPCsByRangeCoord<&Area_t::width>(beginTPCwithArea, endTPCwithArea);
736  std::vector<TPCwithArea_t const*> const& TPCs = TPCgroupInfo.first;
738  TPCgroupInfo.second;
739 
740  if (TPCs.empty()) return {}; // failure?
741  // with only one TPC, then makeTPCPartitionElement() should be used instead!
742  if (TPCs.size() < 4) return {};
743 
744  unsigned int const nWidthParts = TPCgroups.size();
745  if (nWidthParts <= 1) return {}; // only one group ain't no good
746 
747  //
748  // sort TPCs in the first width partition by depth coordinate
749  //
750  auto const FirstColGroupInfo =
751  sortAndGroupTPCsByRangeCoord<&Area_t::depth>(TPCgroups[0], TPCgroups[1]);
752  std::vector<TPCwithArea_t const*> const& FirstColTPCs = FirstColGroupInfo.first;
754  FirstColGroupInfo.second;
755 
756  if (FirstColTPCs.empty()) return {}; // failure?
757  if (FirstColGroups.size() <= 1) return {}; // only one row ain't good either
758 
759  //
760  // collect all candidate separation ranges
761  //
762  // First depth partition has no lower limit, last one has no upper limit
763  // (they include all TPCs with depth lower than the upper limit in the first
764  // case, all TPCs with depth higher of the lower limit in the last case).
765  // Checks need to be done in the gaps between depth partitions.
766  // So we start by skipping the first border.
767 
769  std::vector<Area_t::Range_t> depthGaps; // candidate gaps
770  auto icnext = FirstColGroups.cbegin(), icprev = icnext, icend = FirstColGroups.cend();
771  while (++icnext != icend) {
772  //
773  // establish the range of the group in the candidate depth partition
774  // from the group in the first width partition
775  //
776  auto const cprev = *icprev;
777  auto const cnext = *icnext;
778 
779  depthGaps.emplace_back((*cprev)->area().depth.upper, (*cnext)->area().depth.lower);
780 
781  icprev = icnext;
782  } // while
783  assert(!depthGaps.empty());
784 
785  //
786  // see that for every other width partition separations hold
787  //
788  auto igbegin = TPCgroups.cbegin();
789  while (++igbegin != TPCgroups.cend()) {
790  //
791  // prepare the TPC groups within this width partition
792  //
793  auto igend = std::next(igbegin);
794  auto gbegin = *igbegin;
795  auto gend = (igend == TPCgroups.cend()) ? TPCs.cend() : *igend;
796 
797  auto const ColGroupInfo = sortAndGroupTPCsByRangeCoord<&Area_t::depth>(gbegin, gend);
798  std::vector<TPCwithArea_t const*> const& ColTPCs = ColGroupInfo.first;
800  ColGroupInfo.second;
801 
802  // failure to partition a single column means total failure
803  if (ColTPCs.empty()) return {};
804  if (ColGroups.size() <= 1) return {}; // only one row ain't good either
805 
806  //
807  // compute the coverage of each of the depth groups
808  //
809  std::vector<TPCwithArea_t::Area_t::Range_t> groupDepths(ColGroups.size());
810  auto iGDepth = groupDepths.begin();
811  for (auto icgstart = ColGroups.cbegin(); icgstart != ColGroups.cend(); ++icgstart, ++iGDepth) {
812  auto const icgend = std::next(icgstart);
813  auto ictpc = *icgstart;
814  auto const ictend = (icgend == ColGroups.cend()) ? ColTPCs.cend() : *icgend;
815  while (ictpc != ictend)
816  iGDepth->extendToInclude((*(ictpc++))->area().depth);
817  } // for
818 
819  //
820  // check each of the remaining candidate gaps
821  //
822  auto iGap = depthGaps.begin();
823  while (iGap != depthGaps.end()) {
824  Area_t::Range_t& gap = *iGap;
825 
826  //
827  // check that the gap holds
828  //
829  bool bGoodGap = false;
830  // first TPC starting after the gap (even immediately after):
831  auto iCGroup = std::lower_bound(
832  groupDepths.cbegin(), groupDepths.cend(), gap.upper, SortTPCwithAreaByDepth());
833 
834  // any TPCs before/after this gap?
835  if ((iCGroup != groupDepths.begin()) && (iCGroup != groupDepths.end())) {
836  Area_t::Range_t const& before = *(std::prev(iCGroup));
837  Area_t::Range_t const& after = *iCGroup;
838  Area_t::Range_t const TPCgap{before.upper, after.lower};
839 
840  // correct the gap
841  if (coordIs.strictlySmaller(iGap->lower, TPCgap.lower)) iGap->lower = TPCgap.lower;
842  if (coordIs.strictlyGreater(iGap->upper, TPCgap.upper)) iGap->upper = TPCgap.upper;
843 
844  // if nothing is left, gap is gone
845  bGoodGap = coordIs.nonSmaller(iGap->upper, iGap->lower);
846  } // if TPCs around the gap
847 
848  //
849  // if the gap has been flagged as bad, remove it
850  //
851  if (bGoodGap)
852  ++iGap;
853  else
854  iGap = depthGaps.erase(iGap);
855 
856  } // while (separation)
857 
858  if (depthGaps.empty()) return {}; // no surviving gaps means failure
859 
860  } // while (width partition)
861 
862  //
863  // turn the gaps into separators
864  //
865  std::vector<double> depthSep;
866  std::transform(depthGaps.cbegin(),
867  depthGaps.cend(),
868  std::back_inserter(depthSep),
869  [](auto const& r) { return (r.lower + r.upper) / 2.0; });
870  unsigned int const nDepthParts = depthSep.size() + 1;
871 
872  //
873  // fill the groups with TPCs, and create subpartitions from each of them
874  //
875  geo::part::Partition<geo::TPCGeo const>::Subpartitions_t subparts(nWidthParts * nDepthParts);
876  Area_t totalArea;
877 
878  unsigned int iWidth = 0;
879  for (auto igbegin = TPCgroups.cbegin(); igbegin != TPCgroups.cend(); ++igbegin, ++iWidth) {
880 
881  // sort TPCs in this group (yes, again; this time we don't group just yet)
882  auto igend = std::next(igbegin);
883  auto gbegin = *igbegin;
884  auto gend = (igend == TPCgroups.cend()) ? TPCs.cend() : *igend;
885  std::vector<TPCwithArea_t const*> ColTPCs(gbegin, gend);
886  std::sort(ColTPCs.begin(), ColTPCs.end(), SortTPCareaByAreaRangeLower<&Area_t::depth>());
887 
888  unsigned int iDepth = 0;
889  auto cgstart = ColTPCs.cbegin(), TPCend = ColTPCs.cend();
890  for (double sep : depthSep) {
891 
892  //
893  // collect all TPCs for this partition
894  //
895  // the first TPC that starts *after* the separator:
896  auto cgend = std::upper_bound(cgstart, TPCend, sep, SortTPCwithAreaByDepth());
897  // if we cut out TPCs that were included because of some tolerance,
898  // recover them now
899  while (cgend != cgstart) {
900  auto cglast = std::prev(cgend);
901  if (coordIs.strictlySmaller((*cglast)->area().depth.lower, sep)) break;
902  cgend = cglast;
903  } // while
904  assert(cgstart != cgend); // separator selection should guarantee this
905 
906  //
907  // create and register the partition
908  //
909  auto part = makePartition(cgstart, cgend);
910  if (!part) return {}; // late failure!
911  totalArea.extendToInclude(part->area());
912  subparts[iDepth * nWidthParts + iWidth] = std::move(part);
913 
914  ++iDepth;
915  cgstart = cgend;
916  } // for all depth separators
917 
918  //
919  // collect all the TPCs after the last separator
920  //
921  auto part = makePartition(cgstart, TPCend);
922  if (!part) return {}; // super-late failure!
923  totalArea.extendToInclude(part->area());
924  subparts[iDepth * nWidthParts + iWidth] = std::move(part);
925 
926  } // for all width partitions
927 
928  return std::make_unique<geo::part::GridPartition<geo::TPCGeo const>>(
929  totalArea, std::move(subparts), nWidthParts, nDepthParts);
930 
931 } // makeGridPartition()
TRandom r
Definition: spectrum.C:23
std::unique_ptr< geo::part::Partition< geo::TPCGeo const > > makePartition(BeginIter beginTPCwithArea, EndIter endTPCwithArea)
Provides simple real number checks.
SortTPCareaByAreaRangeLower<&geo::part::AreaOwner::Area_t::depth > SortTPCwithAreaByDepth
intermediate_table::const_iterator const_iterator
std::vector< std::unique_ptr< Partition_t const >> Subpartitions_t
Type of list of subpartitions. It needs to preserve polymorphism.
Definition: Partitions.h:191
TString part[npart]
Definition: Style.C:32
if(nlines<=0)
lar::util::simple_geo::Rectangle< double > Area_t
Type of area covered by the partition.
Definition: Partitions.h:41
template<typename BeginIter , typename EndIter >
std::unique_ptr< geo::part::Partition< geo::TPCGeo const > > makePartition ( BeginIter  beginTPCwithArea,
EndIter  endTPCwithArea 
)

Definition at line 935 of file DriftPartitions.cxx.

References makeDepthPartition(), makeGridPartition(), makeTPCPartitionElement(), makeWidthPartition(), and util::size().

Referenced by geo::buildDriftVolumes(), makeGridPartition(), makePartition(), and makeTPCPartitionElement().

937 {
938  /*
939  * Organizes a list of TPCs in a hierarchical partition.
940  * Three main elements are used:
941  * - single element partition objects: that's the single TPC end point
942  * - TPC groups organised along width
943  * - TPC groups organised along depth
944  *
945  * The procedure is recursively analysing a set of TPCs:
946  * - if the set is actually one TPC only, use a PartitionElement
947  * - attempt partitioning o a grid; if fails:
948  * - attempt partitioning along width:
949  * * determine overlapping groups: a set of TPCs which share part of the
950  * width range
951  * * recurse on each overlapping group with more than one TPC,
952  * attempting a depth partition
953  * * if that fails, bail out since we don't have code to deal with a layout
954  * with areas overlapping on both directions at the same time
955  * * add the single elements and the overlapping groups to the width
956  * partition
957  * - attempt partitioning along height:
958  * * same algorithm as for width
959  * - pick the partitioning with less elements
960  */
961  using value_type = std::remove_reference_t<decltype(*beginTPCwithArea)>;
962  static_assert(std::is_pointer<value_type>() &&
963  std::is_same<std::decay_t<std::remove_pointer_t<value_type>>, TPCwithArea_t>(),
964  "Iterators must point to TPCwithArea_t pointers.");
965 
966  auto const size = std::distance(beginTPCwithArea, endTPCwithArea);
967  if (size == 1) { return makeTPCPartitionElement(**beginTPCwithArea); }
968 
969  auto gPart = makeGridPartition(beginTPCwithArea, endTPCwithArea);
970  if (gPart) return gPart;
971 
972  auto wPart = makeWidthPartition(beginTPCwithArea, endTPCwithArea);
973  auto dPart = makeDepthPartition(beginTPCwithArea, endTPCwithArea);
974 
975  if (wPart) {
976 
977  if (dPart) { // wPart && dPart
978  if (wPart->nParts() < dPart->nParts())
979  return wPart;
980  else
981  return dPart; // slight preference
982  }
983  else { // wPart && !dPart
984  return wPart; // easy choice
985  }
986  }
987  else {
988 
989  if (dPart) { // !wPart && dPart
990  return dPart; // easy choice
991  }
992  else { // !wPart && !dPart
993  return {}; // failure!!
994  }
995  }
996 
997 } // makePartition(Iter)
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:101
std::unique_ptr< geo::part::Partition< geo::TPCGeo const > > makeDepthPartition(BeginIter beginTPCwithArea, EndIter endTPCwithArea)
std::unique_ptr< geo::part::Partition< geo::TPCGeo const > > makeWidthPartition(BeginIter beginTPCwithArea, EndIter endTPCwithArea)
auto makeTPCPartitionElement(TPCwithArea_t const &TPCinfo)
std::unique_ptr< geo::part::Partition< geo::TPCGeo const > > makeGridPartition(BeginIter beginTPCwithArea, EndIter endTPCwithArea)
std::unique_ptr<geo::DriftPartitions::TPCPartition_t> makePartition ( std::vector< TPCwithArea_t > const &  TPCs)

Definition at line 1017 of file DriftPartitions.cxx.

References util::cbegin(), util::cend(), makeCPointerVector(), and makePartition().

1019 {
1020  // TODO use range library instead:
1021  // auto TPCptrs = TPCs | ranges::views::transform(std::addressof);
1022  auto TPCptrs = makeCPointerVector(TPCs);
1023  using std::cbegin;
1024  using std::cend;
1025  return makePartition(cbegin(TPCptrs), cend(TPCptrs));
1026 } // makePartition(coll)
std::unique_ptr< geo::part::Partition< geo::TPCGeo const > > makePartition(BeginIter beginTPCwithArea, EndIter endTPCwithArea)
decltype(auto) constexpr cend(T &&obj)
ADL-aware version of std::cend.
Definition: StdUtils.h:93
auto makeCPointerVector(BeginIter b, EndIter e)
decltype(auto) constexpr cbegin(T &&obj)
ADL-aware version of std::cbegin.
Definition: StdUtils.h:85
template<typename TPCPartitionResultType , geo::part::AreaOwner::AreaRangeMember_t Range, typename BeginIter , typename EndIter , typename SubpartMaker >
std::unique_ptr<geo::part::Partition<geo::TPCGeo const> > makeSortedPartition ( BeginIter  beginTPCwithArea,
EndIter  endTPCwithArea,
SubpartMaker  subpartMaker 
)

Definition at line 633 of file DriftPartitions.cxx.

References computeTotalArea(), and createSubpartitions().

634 {
635  /*
636  * TPCs in input are arranged into a partition split on width direction.
637  * In case of failure, a null pointer is returned.
638  * Do not use this function for a single TPC.
639  * The algorithm is as follows:
640  *
641  * 1. sort the TPCs by width coordinate
642  * 2. for each one, group it with all the following ones overlapping in width
643  * 3. for each group with:
644  * .1 more than one element:
645  * .1 let the group be partitioned on depth
646  * .2 if that fails, we fail
647  * .2 just one element: create a partition element
648  * 4. the resulting partition is the list of one-TPC "groups" and depth-based
649  * subpartitions of many-TPC groups
650  * 5. if the result is a single partition, it must be a depth partition;
651  * then, there is no room for a partition along width, and we declare
652  * failure
653  */
654 
655  //
656  // sort by coordinate and group TPCs; work with pointers for convenience
657  //
658  auto const TPCgroupInfo = sortAndGroupTPCsByRangeCoord<Range>(beginTPCwithArea, endTPCwithArea);
659  std::vector<TPCwithArea_t const*> const& TPCs = TPCgroupInfo.first;
661  TPCgroupInfo.second;
662 
663  if (TPCs.empty()) return {}; // failure?
664 
665  //
666  // for each group, create a subpartition
667  //
668  auto subparts =
669  createSubpartitions(TPCgroups.cbegin(), TPCgroups.cend(), TPCs.cend(), subpartMaker);
670 
671  // if we have grouped everything in a single unit, we have not done any good
672  if (subparts.size() == 1) return {};
673 
674  //
675  // compute the total area (it might have been merged in a previous loop...)
676  //
677  auto totalArea = computeTotalArea(TPCs.cbegin(), TPCs.cend());
678 
679  //
680  // construct and return the final partition
681  //
682  return std::make_unique<TPCPartitionResultType>(totalArea, std::move(subparts));
683 
684 } // makeSortedPartition()
geo::part::Partition< geo::TPCGeo const >::Subpartitions_t createSubpartitions(BeginIter itTPCbegin, EndIter itTPCend, TPCendIter TPCend, SubpartMaker subpartMaker)
geo::part::AreaOwner::Area_t computeTotalArea(BeginIter TPCbegin, EndIter TPCend)
intermediate_table::const_iterator const_iterator
auto makeTPCPartitionElement ( TPCwithArea_t const &  TPCinfo)

Definition at line 602 of file DriftPartitions.cxx.

References geo::part::AreaOwner::area(), makeDepthPartition(), makeGridPartition(), makePartition(), makeWidthPartition(), and TPCwithArea_t::TPC.

Referenced by createSubpartitions(), and makePartition().

603 {
604  return std::make_unique<geo::part::PartitionElement<geo::TPCGeo const>>(TPCinfo.area(),
605  TPCinfo.TPC);
606 }
template<typename BeginIter , typename EndIter >
std::unique_ptr< geo::part::Partition< geo::TPCGeo const > > makeWidthPartition ( BeginIter  beginTPCwithArea,
EndIter  endTPCwithArea 
)

Definition at line 688 of file DriftPartitions.cxx.

References lar::util::simple_geo::Rectangle< double >::width.

Referenced by makePartition(), and makeTPCPartitionElement().

691 {
692  return makeSortedPartition<geo::part::WidthPartition<geo::TPCGeo const>,
694  beginTPCwithArea, endTPCwithArea, &makeDepthPartition<BeginIter, EndIter>);
695 } // makeWidthPartition()
Range_t width
Range along width direction.
Definition: SimpleGeo.h:429
template<geo::part::AreaOwner::AreaRangeMember_t sortingRange, typename BeginIter , typename EndIter >
std::pair<std::vector<TPCwithArea_t const*>, std::vector<std::vector<TPCwithArea_t const*>::const_iterator> > sortAndGroupTPCsByRangeCoord ( BeginIter  beginTPCwithArea,
EndIter  endTPCwithArea 
)

Definition at line 525 of file DriftPartitions.cxx.

526 {
527  //
528  // sort by coordinate; work with pointers for convenience
529  //
530  std::vector<TPCwithArea_t const*> TPCs(beginTPCwithArea, endTPCwithArea);
531  if (TPCs.size() <= 1) return {}; // with only one TPC, refuse to operate
532  std::sort(TPCs.begin(), TPCs.end(), SortTPCareaByAreaRangeLower<sortingRange>());
533 
534  //
535  // group
536  //
538  groupTPCsByRangeCoord<sortingRange>(TPCs.cbegin(), TPCs.cend());
539  assert(!TPCgroups.empty());
540 
541  return {std::move(TPCs), std::move(TPCgroups)};
542 
543 } // sortAndGroupTPCsByRangeCoord()
intermediate_table::const_iterator const_iterator
std::vector<TPCandPos_t> sortTPCsByDriftCoord ( std::vector< geo::TPCGeo const * > const &  TPCs,
geo::DriftPartitions::Decomposer_t const &  decomp 
)

Definition at line 208 of file DriftPartitions.cxx.

References geo::DriftPartitions::driftCoord(), and geo::Decomposer< Vector, Point, ProjVector >::PointNormalComponent().

Referenced by geo::buildDriftVolumes().

210 {
211  /*
212  * TPCs in the argument are sorted by their drift coordinate.
213  * The drift coordinate is defined as the coordinate specified by the normal
214  * direction of the decomposer in argument.
215  * The result is a collection of data structures containing each a TPC and
216  * the drift coordinate that was used to sort it.
217  *
218  * Sorting happens by the drift coordinate of the first wire plane;
219  * the absolute value of the drift coordinate value is not relevant nor it is
220  * well defined.
221  * The result preserves that coordinate for further processing (grouping).
222  */
223  auto const driftCoord = [&decomp](geo::TPCGeo const& TPC) {
224  return decomp.PointNormalComponent(
225  geo::vect::convertTo<geo::DriftPartitions::Position_t>(TPC.FirstPlane().GetCenter()));
226  };
227 
228  std::vector<TPCandPos_t> result;
229  result.reserve(TPCs.size());
230  std::transform(
231  TPCs.cbegin(), TPCs.cend(), std::back_inserter(result), [&driftCoord](geo::TPCGeo const* pTPC) {
232  return TPCandPos_t(driftCoord(*pTPC), pTPC);
233  });
234  // std::pair sorts by first key first, and second key on par
235  // (on par, which may happen often, we don't have means to decide here...)
236  std::sort(result.begin(), result.end());
237  return result;
238 } // sortTPCsByDriftCoord()
Geometry information for a single TPC.
Definition: TPCGeo.h:36
std::pair< double, geo::TPCGeo const * > TPCandPos_t
geo::part::AreaOwner::Area_t TPCarea ( geo::TPCGeo const &  TPC,
geo::DriftPartitions::Decomposer_t const &  decomposer 
)

Definition at line 378 of file DriftPartitions.cxx.

References geo::BoxBoundedGeo::MaxX(), geo::BoxBoundedGeo::MaxY(), geo::BoxBoundedGeo::MaxZ(), geo::BoxBoundedGeo::MinX(), geo::BoxBoundedGeo::MinY(), geo::BoxBoundedGeo::MinZ(), and geo::Decomposer< Vector, Point, ProjVector >::ProjectPointOnPlane().

Referenced by addAreaToTPCs().

380 {
381  /*
382  * Returns the "area" of the TPC.
383  * The area is delimited by TPC bounding box
384  */
385 
386  geo::DriftPartitions::Position_t const lower{TPC.MinX(), TPC.MinY(), TPC.MinZ()};
387  geo::DriftPartitions::Position_t const upper{TPC.MaxX(), TPC.MaxY(), TPC.MaxZ()};
388 
389  auto const lowerProj = decomposer.ProjectPointOnPlane(lower);
390  auto const upperProj = decomposer.ProjectPointOnPlane(upper);
391 
392  // we ask to sort the ranges, since the reference base may be flipped
393  return {{lowerProj.X(), upperProj.X(), true}, {lowerProj.Y(), upperProj.Y(), true}};
394 
395 } // TPCarea()
geo::Point_t Position_t
Type representing a position in 3D space.