src/segmentations/Segmentation.cpp

0001 //==========================================================================
0002 //  AIDA Detector description implementation
0003 //--------------------------------------------------------------------------
0004 // Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
0005 // All rights reserved.
0006 //
0007 // For the licensing terms see $DD4hepINSTALL/LICENSE.
0008 // For the list of contributors see $DD4hepINSTALL/doc/CREDITS.
0009 //
0010 //==========================================================================
0011 /*
0012  * Segmentation.cpp
0013  *
0014  *  Created on: Jun 27, 2013
0015  *      Author: Christian Grefe, CERN
0016  */
0017
0018 #include <DDSegmentation/Segmentation.h>
0019
0020 #include <iostream>
0021 #include <sstream>
0022 #include <stdexcept>
0023 #include <cmath>
0024 #include <algorithm>
0025 #include <iomanip>
0026
0027 namespace dd4hep {
0028
0029   namespace DDSegmentation {
0030
0031     using std::cerr;
0032     using std::endl;
0033     using std::map;
0034     using std::runtime_error;
0035     using std::stringstream;
0036     using std::vector;
0037
0038     /// Default constructor used by derived classes passing the encoding string
0039     Segmentation::Segmentation(const std::string& cellEncoding) :
0040       _name("Segmentation"), _type("Segmentation"), _decoder(new BitFieldCoder(cellEncoding)), _ownsDecoder(true) {
0041
0042     }
0043
0044     /// Default constructor used by derived classes passing an existing decoder
0045     Segmentation::Segmentation(const BitFieldCoder* newDecoder) :
0046       _name("Segmentation"), _type("Segmentation"), _decoder(newDecoder), _ownsDecoder(false) {
0047     }
0048
0049     /// Destructor
0050     Segmentation::~Segmentation() {
0051       if (_ownsDecoder and _decoder != 0) {
0052         delete _decoder;
0053       }
0054       for (auto& p : _parameters)  {
0055         if ( p.second ) delete p.second;
0056       }
0057       _parameters.clear();
0058     }
0059
0060     /// Add subsegmentation. Call only valid for Multi-segmentations. Default implementation throws an exception
0061     void Segmentation::addSubsegmentation(long /* key_min */, long /* key_max */, Segmentation* /* entry */)   {
0062       throw std::runtime_error("This segmentation type:"+_type+" does not support sub-segmentations.");
0063     }
0064
0065     /// Determine the volume ID from the full cell ID by removing all local fields
0066     VolumeID Segmentation::volumeID(const CellID& cID) const {
0067       map<std::string, StringParameter>::const_iterator it;
0068       VolumeID vID = cID ;
0069       for (it = _indexIdentifiers.begin(); it != _indexIdentifiers.end(); ++it) {
0070         std::string identifier = it->second->typedValue();
0071         _decoder->set(vID,identifier,0);
0072       }
0073       return vID;
0074     }
0075
0076     /// Calculates the neighbours of the given cell ID and adds them to the list of neighbours
0077     void Segmentation::neighbours(const CellID& cID, std::set<CellID>& cellNeighbours) const {
0078       map<std::string, StringParameter>::const_iterator it;
0079       for (it = _indexIdentifiers.begin(); it != _indexIdentifiers.end(); ++it) {
0080         const std::string& identifier = it->second->typedValue();
0081         CellID nID = cID ;
0082         int currentValue = _decoder->get(cID,identifier);
0083         // add both neighbouring cell IDs, don't add out of bound indices
0084         try {
0085           _decoder->set(nID,identifier,currentValue - 1);
0086           cellNeighbours.insert(nID);
0087         } catch (runtime_error& e) {
0088           // nothing to do
0089         }
0090         try {
0091           _decoder->set(nID,identifier,currentValue + 1);
0092           cellNeighbours.insert(nID);
0093         } catch (runtime_error& e) {
0094           // nothing to do
0095         }
0096       }
0097     }
0098
0099     /// Set the underlying decoder
0100     void Segmentation::setDecoder(const BitFieldCoder* newDecoder) {
0101       if ( _decoder == newDecoder )
0102         return; //self assignment
0103       else if (_ownsDecoder)
0104         delete _decoder;
0105       _decoder = newDecoder;
0106       _ownsDecoder = false;
0107     }
0108
0109     /// Access to parameter by name
0110     Parameter Segmentation::parameter(const std::string& parameterName) const {
0111       map<std::string, Parameter>::const_iterator it = _parameters.find(parameterName);
0112       if (it != _parameters.end()) {
0113         return it->second;
0114       }
0115       stringstream s;
0116       s << "Unknown parameter " << parameterName << " for segmentation type " << _type;
0117       throw std::runtime_error(s.str());
0118     }
0119
0120     /// Access to all parameters
0121     Parameters Segmentation::parameters() const {
0122       Parameters pars;
0123       for ( const auto& it : _parameters )
0124         pars.emplace_back(it.second);
0125       return pars;
0126     }
0127
0128     /// Set all parameters from an existing set of parameters
0129     void Segmentation::setParameters(const Parameters& pars) {
0130       for ( const auto* p : pars )
0131         parameter(p->name())->value() = p->value();
0132     }
0133
0134     /// Add a cell identifier to this segmentation. Used by derived classes to define their required identifiers
0135     void Segmentation::registerIdentifier(const std::string& idName, const std::string& idDescription, std::string& identifier,
0136                                           const std::string& defaultValue) {
0137       StringParameter idParameter =
0138         new TypedSegmentationParameter<std::string>(idName, idDescription, identifier, defaultValue,
0139                                                     SegmentationParameter::NoUnit, true);
0140       _parameters[idName]       = idParameter;
0141       _indexIdentifiers[idName] = idParameter;
0142     }
0143
0144     /// Helper method to convert a bin number to a 1D position
0145     double Segmentation::binToPosition(FieldID bin, double cellSize, double offset) {
0146       return bin * cellSize + offset;
0147     }
0148
0149     /// Helper method to convert a 1D position to a cell ID
0150     int Segmentation::positionToBin(double position, double cellSize, double offset) {
0151       if (cellSize <= 1e-10) {
0152         throw runtime_error("Invalid cell size: 0.0");
0153       }
0154       return int(floor((position + 0.5 * cellSize - offset) / cellSize));
0155     }
0156
0157     /// Helper method to convert a bin number to a 1D position given a vector of binBoundaries
0158     double Segmentation::binToPosition(FieldID bin, std::vector<double> const& cellBoundaries, double offset) {
0159       return (cellBoundaries[bin+1] + cellBoundaries[bin])*0.5 + offset;
0160     }
0161     /// Helper method to convert a 1D position to a cell ID given a vector of binBoundaries
0162     int Segmentation::positionToBin(double position, std::vector<double> const& cellBoundaries, double offset) {
0163
0164       // include the lower edge to the segmentation, deal with numerical issues
0165       if(fabs(position/cellBoundaries.front()-1.0) < 3e-12) return 0;
0166
0167       // include the upper edge of the last bin to the segmentation, deal with numerical issues
0168       if(fabs(position/cellBoundaries.back()-1.0)  < 3e-12) return int(cellBoundaries.size()-2);
0169
0170       // hits outside cannot be treated
0171       if(position < cellBoundaries.front()) {
0172         std::stringstream err;
0173         err << std::setprecision(20) << std::scientific;
0174         err << "Hit Position (" << position << ") is below the acceptance"
0175             << " (" << cellBoundaries.front() << ") "
0176             << "of the segmentation";
0177         throw std::runtime_error(err.str());
0178       }
0179       if(position > cellBoundaries.back() ) {
0180         std::stringstream err;
0181         err << std::setprecision(20) << std::scientific;
0182         err << "Hit Position (" << position << ") is above the acceptance"
0183             << " (" << cellBoundaries.back() << ") "
0184             << "of the segmentation";
0185         throw std::runtime_error(err.str());
0186       }
0187       std::vector<double>::const_iterator bin = std::upper_bound(cellBoundaries.begin(),
0188                                                                  cellBoundaries.end(),
0189                                                                  position-offset);
0190       // need to reduce found bin by one, because upper_bound works that way, lower_bound would not help
0191       return bin - cellBoundaries.begin() - 1 ;
0192
0193     }
0194
0195     std::vector<double> Segmentation::cellDimensions(const CellID&) const {
0196       std::stringstream errorMessage;
0197       errorMessage << __func__ << " is not implemented for " << _type;
0198       throw std::logic_error(errorMessage.str());
0199     }
0200
0201   } /* namespace DDSegmentation */
0202 } /* namespace dd4hep */