diff --git a/DDRec/include/DDRec/API/Exceptions.h b/DDRec/include/DDRec/API/Exceptions.h index 8ad8e12a870b601ca29671561e309cb8ea68394b..24a76f78a4e31b310f8ac6aea4b2ed74ce009c73 100644 --- a/DDRec/include/DDRec/API/Exceptions.h +++ b/DDRec/include/DDRec/API/Exceptions.h @@ -30,7 +30,7 @@ public: private: static std::string createMsg(const std::string& msg, const DDSegmentation::CellID& cellID) { std::stringstream s; - s << msg; + s << "Invalid cell ID: " << msg; if (cellID) s << " (" << cellID << ")"; return s.str(); @@ -45,7 +45,7 @@ public: private: static std::string createMsg(const std::string& msg, const Geometry::Position& position) { std::stringstream s; - s << msg << " " << position; + s << "Invalid position: " << msg << " " << position; return s.str(); } }; @@ -58,7 +58,7 @@ public: private: static std::string createMsg(const std::string& msg, const Geometry::DetElement& det) { std::stringstream s; - s << msg; + s << "Invalid detector element: " << msg; if (det.isValid()) s << " (" <<det.name() << ")"; return s.str(); diff --git a/DDRec/include/DDRec/API/IDDecoder.h b/DDRec/include/DDRec/API/IDDecoder.h index 1d555b1d99b6567c8a833eae80d0e297c8856038..96358f5ad78cd6387f2bff8b08eeeff6e8e551e2 100644 --- a/DDRec/include/DDRec/API/IDDecoder.h +++ b/DDRec/include/DDRec/API/IDDecoder.h @@ -31,98 +31,137 @@ typedef DDSegmentation::VolumeID VolumeID; */ class IDDecoder { public: - /** - * Default constructor using the name of the corresponding readout collection - */ - IDDecoder(const std::string& collectionName); - - /** - * Default constructor using a readout object - */ - IDDecoder(const Geometry::Readout& readout); - - /** - * Destructor - */ - virtual ~IDDecoder(); - - /** - * Returns the cell ID from the local position in the given volume ID. - */ + class BarrelEndcapFlag { + public: + enum BarrelEncapID { + Barrel=0, + EndcapSouth, + EndcapNorth, + n_BarrelEndcapID + }; + BarrelEndcapFlag(unsigned int value) : + value(static_cast<BarrelEncapID>(value)) {} + + virtual ~BarrelEndcapFlag() {} + + bool isBarrel() const { + return value == Barrel; + } + + bool isEndcap() const { + return value == EndcapNorth || value == EndcapSouth; + } + + bool isEndcapSouth() const { + return value == EndcapSouth; + } + + bool isEndcapNorth() const { + return value == EndcapSouth; + } + + BarrelEncapID getValue() const { + return value; + } + + protected: + BarrelEncapID value; + }; + + /// Access to the global IDDecoder instance + static IDDecoder& getInstance(); + + /// Destructor + virtual ~IDDecoder() {}; + + /// Returns the cell ID from the local position in the given volume ID. CellID cellIDFromLocal(const Geometry::Position& local, const VolumeID volumeID) const; - /** - * Returns the global cell ID from a given global position - */ + /// Returns the global cell ID from a given global position CellID cellID(const Geometry::Position& global) const; - /** - * Returns the global position from a given cell ID - */ + /// Returns the global position from a given cell ID Geometry::Position position(const CellID& cellID) const; - /* - * Returns the local position from a given cell ID - */ + /// Returns the local position from a given cell ID Geometry::Position localPosition(const CellID& cellID) const; - /* - * Returns the volume ID of a given cell ID - */ + /// Returns the volume ID of a given cell ID VolumeID volumeID(const CellID& cellID) const; - /* - * Returns the volume ID of a given global position - */ + /// Returns the volume ID of a given global position VolumeID volumeID(const Geometry::Position& global) const; - /* - * Returns the placement for a given cell ID - */ + /// Returns the placement for a given cell ID Geometry::PlacedVolume placement(const CellID& cellID) const; - /* - * Returns the placement for a given global position - */ + /// Returns the placement for a given global position Geometry::PlacedVolume placement(const Geometry::Position& global) const; - /* - * Returns the subdetector for a given cell ID - */ + /// Returns the subdetector for a given cell ID Geometry::DetElement subDetector(const CellID& cellID) const; - /* - * Returns the subdetector for a given global position - */ + /// Returns the subdetector for a given global position Geometry::DetElement subDetector(const Geometry::Position& global) const; - /* - * Returns the closest detector element in the hierarchy for a given cell ID - */ + /// Returns the closest detector element in the hierarchy for a given cell ID Geometry::DetElement detectorElement(const CellID& cellID) const; - /* - * Returns the closest detector element in the hierarchy for a given global position - */ + /// Returns the closest detector element in the hierarchy for a given global position Geometry::DetElement detectorElement(const Geometry::Position& global) const; - /* - * Calculates the neighbours of the given cell ID and adds them to the list of neighbours - */ + /// Access to the Readout object for a given cell ID + Geometry::Readout readout(const CellID& cellID) const; + + /// Access to the Readout object for a given global position + Geometry::Readout readout(const Geometry::Position& global) const; + + /// Calculates the neighbours of the given cell ID and adds them to the list of neighbours void neighbours(const CellID& cellID, std::set<CellID>& neighbours) const; - /* - * Checks if the given cell IDs are neighbours - */ + /// Checks if the given cell IDs are neighbours bool areNeighbours(const CellID& cellID, const CellID& otherCellID) const; + /// Access to the barrel-endcap flag + BarrelEndcapFlag barrelEndcapFlag(const CellID& cellID) const; + + /// Access to the layer index + long int layerIndex(const CellID& cellID) const; + + /// Access to the system index + long int systemIndex(const CellID& cellID) const; + + static std::string barrelIdentifier() { + return std::string("barrel"); + } + + static std::string layerIdentifier() { + return std::string("layer"); + } + + static std::string systemIdentifier() { + return std::string("system"); + } + protected: - Geometry::Readout _readout; Geometry::VolumeManager _volumeManager; TGeoManager* _tgeoMgr; - // helper method to get the closest daughter DetElement to the position starting from the given DetElement + /// Helper method to find the corresponding Readout object to a DetElement + Geometry::Readout findReadout(const Geometry::DetElement& det) const; + + /// Helper method to get the closest daughter DetElement to the position starting from the given DetElement static Geometry::DetElement getClosestDaughter(const Geometry::DetElement& det, const Geometry::Position& position); + +private: + /// Default constructor + IDDecoder(); + + /// Disable copy constructor + IDDecoder(const IDDecoder&); + + /// Disable assignment operator + void operator=(const IDDecoder&); }; } /* namespace DDRec */ diff --git a/DDRec/include/DDRec/Extensions/SubdetectorExtensionImpl.h b/DDRec/include/DDRec/Extensions/SubdetectorExtensionImpl.h index 0a37fa2b33d108e400820378683c5fa176e7acb2..3de83e50dc92280cba4c56778966d03c674b0627 100644 --- a/DDRec/include/DDRec/Extensions/SubdetectorExtensionImpl.h +++ b/DDRec/include/DDRec/Extensions/SubdetectorExtensionImpl.h @@ -84,7 +84,9 @@ public: protected: Geometry::DetElement det; bool _isBarrel; + bool _setIsBarrel; bool _isEndcap; + bool _setIsEndcap; double _rMin; bool _setRMin; double _rMax; diff --git a/DDRec/src/IDDecoder.cpp b/DDRec/src/IDDecoder.cpp index 9c8b1ce48f08ff433a6b2ce46a57fce93379fbe7..213eb374f502f3499fc334128a6ce7d524dbb1e2 100644 --- a/DDRec/src/IDDecoder.cpp +++ b/DDRec/src/IDDecoder.cpp @@ -6,6 +6,8 @@ */ #include "DDRec/API/IDDecoder.h" +#include "DDRec/API/Exceptions.h" + #include "DD4hep/LCDD.h" #include "DD4hep/VolumeManager.h" @@ -22,36 +24,19 @@ using Geometry::VolumeManager; using Geometry::Volume; using std::set; -/** - * Default constructor using the name of the corresponding readout collection - */ -IDDecoder::IDDecoder(const std::string& collectionName) { - LCDD& lcdd = LCDD::getInstance(); - _readout = lcdd.readout(collectionName); - _volumeManager = lcdd.volumeManager(); - if (not _volumeManager.isValid()) { - _volumeManager = VolumeManager(lcdd, "volman", lcdd.world(), Readout(), VolumeManager::TREE); - } - _tgeoMgr = Geometry::LCDD::getInstance().world().volume()->GetGeoManager(); +IDDecoder& IDDecoder::getInstance() { + static IDDecoder idd; + return idd; } -/** - * Default constructor using a readout object - */ -IDDecoder::IDDecoder(const Readout& readout) { +/// Default constructor +IDDecoder::IDDecoder() { LCDD& lcdd = LCDD::getInstance(); - _readout = readout; _volumeManager = lcdd.volumeManager(); if (not _volumeManager.isValid()) { _volumeManager = VolumeManager(lcdd, "volman", lcdd.world(), Readout(), VolumeManager::TREE); } - _tgeoMgr = Geometry::LCDD::getInstance().world().volume()->GetGeoManager(); -} - -/** - * Destructor - */ -IDDecoder::~IDDecoder() { + _tgeoMgr = lcdd.world().volume()->GetGeoManager(); } /** @@ -63,10 +48,11 @@ CellID IDDecoder::cellIDFromLocal(const Position& local, const VolumeID volumeID local.GetCoordinates(l); // FIXME: direct lookup of transformations seems to be broken //const TGeoMatrix& localToGlobal = _volumeManager.worldTransformation(volumeID); - const TGeoMatrix& localToGlobal = this->detectorElement(volumeID).worldTransformation(); + DetElement det = this->detectorElement(volumeID); + const TGeoMatrix& localToGlobal = det.worldTransformation(); localToGlobal.LocalToMaster(l, g); Position global(g[0], g[1], g[2]); - return _readout.segmentation().cellID(local, global, volumeID); + return this->findReadout(det).segmentation().cellID(local, global, volumeID); } /** @@ -79,10 +65,11 @@ CellID IDDecoder::cellID(const Position& global) const { global.GetCoordinates(g); // FIXME: direct lookup of transformations seems to be broken //const TGeoMatrix& localToGlobal = _volumeManager.worldTransformation(volID); - const TGeoMatrix& localToGlobal = this->detectorElement(volID).worldTransformation(); + DetElement det = this->detectorElement(volID); + const TGeoMatrix& localToGlobal = det.worldTransformation(); localToGlobal.MasterToLocal(g, l); Position local(l[0], l[1], l[2]); - return _readout.segmentation().cellID(local, global, volID); + return this->findReadout(det).segmentation().cellID(local, global, volID); } /** @@ -91,11 +78,12 @@ CellID IDDecoder::cellID(const Position& global) const { Position IDDecoder::position(const CellID& cellID) const { double l[3]; double g[3]; - Position local = _readout.segmentation().position(cellID); + DetElement det = this->detectorElement(cellID); + Position local = this->findReadout(det).segmentation().position(cellID); local.GetCoordinates(l); // FIXME: direct lookup of transformations seems to be broken //const TGeoMatrix& localToGlobal = _volumeManager.worldTransformation(cellID); - const TGeoMatrix& localToGlobal = this->detectorElement(cellID).worldTransformation(); + const TGeoMatrix& localToGlobal = det.worldTransformation(); localToGlobal.LocalToMaster(l, g); return Position(g[0], g[1], g[2]); } @@ -104,21 +92,23 @@ Position IDDecoder::position(const CellID& cellID) const { * Returns the local position from a given cell ID */ Position IDDecoder::localPosition(const CellID& cellID) const { - return _readout.segmentation().position(cellID); + DetElement det = this->detectorElement(cellID); + return this->findReadout(det).segmentation().position(cellID); } /* * Returns the volume ID of a given cell ID */ VolumeID IDDecoder::volumeID(const CellID& cellID) const { - return _readout.segmentation()->volumeID(cellID); + DetElement det = this->detectorElement(cellID); + return this->findReadout(det).segmentation()->volumeID(cellID); } /* * Returns the volume ID of a given global position */ VolumeID IDDecoder::volumeID(const Position& position) const { - DetElement det = detectorElement(position); + DetElement det = this->detectorElement(position); return det.volumeID(); } @@ -161,17 +151,33 @@ DetElement IDDecoder::detectorElement(const CellID& cellID) const { * Returns the closest detector element in the hierarchy for a given global position */ DetElement IDDecoder::detectorElement(const Position& position) const { - Geometry::DetElement world = Geometry::LCDD::getInstance().world(); - Geometry::DetElement det = getClosestDaughter(world, position); - // Fixme: check if this is valid, otherwise throw exception + DetElement world = Geometry::LCDD::getInstance().world(); + DetElement det = getClosestDaughter(world, position); + if (not det.isValid()) { + throw invalid_position("DD4hep::DDRec::IDDecoder::detectorElement", position); + } + std::cout << det.name() << std::endl; return det; } +/// Access to the Readout object for a given cell ID +Geometry::Readout IDDecoder::readout(const CellID& cellID) const { + DetElement det = this->detectorElement(cellID); + return this->findReadout(det); +} + +/// Access to the Readout object for a given global position +Geometry::Readout IDDecoder::readout(const Geometry::Position& global) const { + DetElement det = this->detectorElement(global); + return this->findReadout(det); +} + /* * Calculates the neighbours of the given cell ID and adds them to the list of neighbours */ void IDDecoder::neighbours(const CellID& cellID, set<CellID>& neighbours) const { - _readout.segmentation()->neighbours(cellID, neighbours); + DetElement det = this->detectorElement(cellID); + this->findReadout(det).segmentation()->neighbours(cellID, neighbours); } /* @@ -179,13 +185,58 @@ void IDDecoder::neighbours(const CellID& cellID, set<CellID>& neighbours) const */ bool IDDecoder::areNeighbours(const CellID& cellID, const CellID& otherCellID) const { set<CellID> neighbours; - _readout.segmentation()->neighbours(cellID, neighbours); + DetElement det = this->detectorElement(cellID); + this->findReadout(det).segmentation()->neighbours(cellID, neighbours); return neighbours.count(otherCellID) != 0; } -Geometry::DetElement IDDecoder::getClosestDaughter(const Geometry::DetElement& det, - const Geometry::Position& position) { - Geometry::DetElement result; +/// Access to the barrel-endcap flag +IDDecoder::BarrelEndcapFlag IDDecoder::barrelEndcapFlag(const CellID& cellID) const { + Readout r = this->readout(cellID); + return BarrelEndcapFlag(r.idSpec().field(this->barrelIdentifier())->value(cellID)); +} + +/// Access to the layer index +long int IDDecoder::layerIndex(const CellID& cellID) const { + Readout r = this->readout(cellID); + return r.idSpec().field(this->layerIdentifier())->value(cellID); +} + +/// Access to the system index +long int IDDecoder::systemIndex(const CellID& cellID) const { + Readout r = this->readout(cellID); + return r.idSpec().field(this->systemIdentifier())->value(cellID); +} + +// helper method to find the corresponding Readout object to a DetElement +Readout IDDecoder::findReadout(const Geometry::DetElement& det) const { + + // first check if top level is a sensitive detector + if (det.volume().isValid() and det.volume().isSensitive()) { + Geometry::SensitiveDetector sd = det.volume().sensitiveDetector(); + if (sd.isValid() and sd.readout().isValid()) { + return sd.readout(); + } + } + + // check all children recursively for the first valid Readout object + const DetElement::Children& children = det.children(); + DetElement::Children::const_iterator it = children.begin(); + while (it != children.end()) { + Readout r = findReadout(it->second); + if (r.isValid()) { + return r; + } + ++it; + } + + // neither this or any daughter is sensitive + return Readout(); +} + +// helper method to get the closest daughter DetElement to the position starting from the given DetElement +DetElement IDDecoder::getClosestDaughter(const DetElement& det, const Position& position) { + DetElement result; // check if we have a shape and see if we are inside if (det.volume().isValid() and det.volume().solid().isValid()) { @@ -196,14 +247,14 @@ Geometry::DetElement IDDecoder::getClosestDaughter(const Geometry::DetElement& d result = det; } else { // assuming that any daughter shape would be inside this shape - return Geometry::DetElement(); + return DetElement(); } } const DetElement::Children& children = det.children(); DetElement::Children::const_iterator it = children.begin(); while (it != children.end()) { - Geometry::DetElement daughterDet = getClosestDaughter(it->second, position); + DetElement daughterDet = getClosestDaughter(it->second, position); if (daughterDet.isValid()) { result = daughterDet; break; diff --git a/DDRec/src/SubdetectorExtensionImpl.cpp b/DDRec/src/SubdetectorExtensionImpl.cpp index f74bd1c4f97a2b79fdf4707d26a780eb25dfbc3b..296d05ae4605f3b51df33f6b491d2de0c84299b4 100644 --- a/DDRec/src/SubdetectorExtensionImpl.cpp +++ b/DDRec/src/SubdetectorExtensionImpl.cpp @@ -6,6 +6,7 @@ */ #include "DDRec/Extensions/SubdetectorExtensionImpl.h" +#include "DDRec/API/IDDecoder.h" #include "DD4hep/LCDD.h" @@ -23,8 +24,12 @@ SubdetectorExtensionImpl::SubdetectorExtensionImpl(const Geometry::DetElement& d /// Copy constructor SubdetectorExtensionImpl::SubdetectorExtensionImpl(const SubdetectorExtensionImpl& e, const Geometry::DetElement& d) { this->det = d; - this->setIsBarrel(e.isBarrel()); - this->setIsEndcap(e.isEndcap()); + if (e._setIsBarrel) { + this->setIsBarrel(e.isBarrel()); + } + if (e._setIsEndcap) { + this->setIsEndcap(e.isEndcap()); + } if (e._setRMin) { this->setRMin(e.getRMin()); } @@ -49,14 +54,20 @@ SubdetectorExtensionImpl::~SubdetectorExtensionImpl() { /// Is this a barrel detector bool SubdetectorExtensionImpl::isBarrel() const { - /// Fixme: could use IDDecoder to figure this out - return _isBarrel; + if (_setIsBarrel) { + return _isBarrel; + } + IDDecoder& id = IDDecoder::getInstance(); + return id.barrelEndcapFlag(det.volumeID()).isBarrel(); } /// Is this an endcap detector bool SubdetectorExtensionImpl::isEndcap() const { - /// Fixme: could use IDDecoder to figure this out - return _isEndcap; + if (_setIsEndcap) { + return _isEndcap; + } + IDDecoder& id = IDDecoder::getInstance(); + return id.barrelEndcapFlag(det.volumeID()).isEndcap(); } /// Access to the inner radius @@ -66,9 +77,19 @@ double SubdetectorExtensionImpl::getRMin() const { } if (det.isValid() and det.volume().isValid() and det.volume().solid().isValid()) { Geometry::Solid solid = det.volume().solid(); - Geometry::Tube tube(solid); - if (tube.isValid()) { - return tube->GetRmin(); + try { + Geometry::Tube tube(solid); + if (tube.isValid()) { + return tube->GetRmin(); + } + } catch (std::runtime_error& e) { + } + try { + Geometry::PolyhedraRegular polyhedra(solid); + if (polyhedra.isValid()) { + return polyhedra->GetRmin()[0]; + } + } catch (std::runtime_error& e) { } } return 0.; @@ -81,9 +102,19 @@ double SubdetectorExtensionImpl::getRMax() const { } if (det.isValid() and det.volume().isValid() and det.volume().solid().isValid()) { Geometry::Solid solid = det.volume().solid(); - Geometry::Tube tube(solid); - if (tube.isValid()) { - return tube->GetRmax(); + try { + Geometry::Tube tube(solid); + if (tube.isValid()) { + return tube->GetRmin(); + } + } catch (std::runtime_error& e) { + } + try { + Geometry::PolyhedraRegular polyhedra(solid); + if (polyhedra.isValid()) { + return polyhedra->GetRmax()[0]; + } + } catch (std::runtime_error& e) { } } return 0.; @@ -98,7 +129,10 @@ double SubdetectorExtensionImpl::getZMin() const { Geometry::Solid solid = det.volume().solid(); Geometry::Box box(solid); if (box.isValid()) { - return box->GetDZ(); + Geometry::Position local(0.,0.,-box->GetDZ()/2.); + Geometry::Position global; + det.localToWorld(local, global); + return global.z(); } } return 0.; @@ -113,7 +147,10 @@ double SubdetectorExtensionImpl::getZMax() const { Geometry::Solid solid = det.volume().solid(); Geometry::Box box(solid); if (box.isValid()) { - return box->GetDZ(); + Geometry::Position local(0.,0.,box->GetDZ()/2.); + Geometry::Position global; + det.localToWorld(local, global); + return global.z(); } } return 0.; @@ -129,9 +166,13 @@ int SubdetectorExtensionImpl::getNSides() const { } if (det.isValid() and det.volume().isValid() and det.volume().solid().isValid()) { Geometry::Solid solid = det.volume().solid(); - Geometry::PolyhedraRegular polyhedra(solid); - if (polyhedra.isValid()) { - return polyhedra->GetNedges(); + try { + Geometry::PolyhedraRegular polyhedra(solid); + if (polyhedra.isValid()) { + return polyhedra->GetNedges(); + } + } catch (std::runtime_error& e) { + } } return 0; @@ -184,7 +225,9 @@ void SubdetectorExtensionImpl::setNSides(int value) { void SubdetectorExtensionImpl::resetAll() { _isBarrel = false; + _setIsBarrel = false; _isEndcap = false; + _setIsEndcap = false; _rMin = 0.; _setRMin = false; _rMax = 0.;