Volumes.cpp

Go to the documentation of this file.

//==========================================================================
//  AIDA Detector description implementation 
//--------------------------------------------------------------------------
// Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
// All rights reserved.
//
// For the licensing terms see $DD4hepINSTALL/LICENSE.
// For the list of contributors see $DD4hepINSTALL/doc/CREDITS.
//
// Author     : M.Frank
//
//==========================================================================
 
// Framework include files
#include <DD4hep/Detector.h>
#include <DD4hep/Printout.h>
#include <DD4hep/InstanceCount.h>
#include <DD4hep/MatrixHelpers.h>
#include <DD4hep/detail/ObjectsInterna.h>
 
// ROOT include files
#include <TROOT.h>
#include <TClass.h>
#include <TColor.h>
#include <TGeoShape.h>
#include <TGeoVolume.h>
#include <TGeoNode.h>
#include <TGeoMatrix.h>
#include <TGeoMedium.h>
 
#include <TGeoVoxelFinder.h>
#include <TGeoShapeAssembly.h>
#include <TGeoScaledShape.h>
#include <TMap.h>
 
// C/C++ include files
#include <climits>
#include <iostream>
#include <stdexcept>
#include <sstream>
#include <iomanip>
 
using namespace dd4hep;
 
/*
 *  The section below uses the new ROOT features using user extensions to volumes
 *  and placements. Once this is common, the above mechanism should be thrown away....
 *
 *  M.Frank
 */
 
typedef TGeoNode                geo_node_t;
typedef TGeoVolume              geo_volume_t;
typedef TGeoVolumeAssembly      geo_assembly_t;
 
ClassImp(VolumeExtension)
ClassImp(PlacedVolumeExtension)
ClassImp(PlacedVolumeExtension::Parameterisation)
 
namespace {
 
  static constexpr double s_rotation_test_limit = 1e-12;
  static bool s_verifyCopyNumbers = true;
 
  template <typename T> typename T::Object* _userExtension(const T& v)  {
    typedef typename T::Object O;
    O* o = (O*)(v.ptr()->GetUserExtension());
    return o;
  }
 
  TGeoVolume* _createTGeoVolume(const std::string& name, TGeoShape* s, TGeoMedium* m)  {
    geo_volume_t* e = new geo_volume_t(name.c_str(),s,m);
    e->SetUserExtension(new Volume::Object());
    return e;
  }
  TGeoVolume* _createTGeoVolumeAssembly(const std::string& name)  {
    geo_assembly_t* e = new geo_assembly_t(name.c_str()); // It is important to use the correct constructor!!
    e->SetUserExtension(new Assembly::Object());
    return e;
  }
  TGeoVolumeMulti* _createTGeoVolumeMulti(const std::string& name, TGeoMedium* medium)  {
    TGeoVolumeMulti* e = new TGeoVolumeMulti(name.c_str(), medium);
    e->SetUserExtension(new VolumeMulti::Object());
    return e;
  }
  PlacedVolume::Object* _data(const PlacedVolume& v) {
    PlacedVolume::Object* o = _userExtension(v);
    if (o) return o;
    except("PlacedVolume::_data", "+++ Attempt to access invalid handle of type: PlacedVolume");
    return nullptr;
  }
  Volume::Object* _data(const Volume& v, bool throw_exception = true) {
    //if ( v.ptr() && v.ptr()->IsA() == TGeoVolume::Class() ) return v.data<Volume::Object>();
    Volume::Object* o = _userExtension(v);
    if (o)
      return o;
    else if (!throw_exception)
      return nullptr;
    except("Volume::_data", "+++ Attempt to access invalid handle of type: Volume");
    return nullptr;
  }
 
  class VolumeImport   {
  public:
    void setShapeTitle(TGeoVolume* vol)   {
      if ( vol )   {
        TGeoShape* sh = vol->GetShape();
        std::string tag = get_shape_tag(sh);
        sh->SetTitle(tag.c_str());
      }
    }
 
  public:
    void operator()(TGeoVolume* v)   {
      TClass* c = v->IsA();
      if ( !v->GetUserExtension() )   {
        if ( c == geo_volume_t::Class() )
          v->SetUserExtension(new Volume::Object());
        else if  ( c == geo_assembly_t::Class() )
          v->SetUserExtension(new Assembly::Object());
        else if  ( c == TGeoVolumeMulti::Class() )
          v->SetUserExtension(new VolumeMulti::Object());
        else
          except("dd4hep","VolumeImport: Unknown TGeoVolume sub-class: %s",v->IsA()->GetName());
      }
      if  ( c == TGeoVolumeMulti::Class() )  {
        TGeoVolumeMulti* mv = (TGeoVolumeMulti*)v;
        for( int i=0, n=mv->GetNvolumes(); i<n; ++i )   {
          TGeoVolume* vol = mv->GetVolume(i);
          this->setShapeTitle(vol);
          (*this)(vol);
        }
      }
      for( Int_t i=0; i<v->GetNdaughters(); ++i )  {
        geo_node_t* pv = v->GetNode(i);
        if ( !pv->GetUserExtension() )
          pv->geo_node_t::SetUserExtension(new PlacedVolume::Object());
        (*this)(pv->GetVolume());
      }
    }
    void operator()(TGeoVolume* new_v, TGeoVolume* old_v, SensitiveDetector sd, int set_bit)   {
      if ( !new_v || !old_v )  {
        except("dd4hep","VolumeImport: ERROR: The refected volume is INVALID!");        
      }
      else if ( !old_v->GetUserExtension() )   {
        except("dd4hep","VolumeImport: ERROR: Reflection of non-dd4hep volume %s",new_v->IsA()->GetName());
      }
      else if ( !new_v->GetUserExtension() )   {
        TClass* c = new_v->IsA();
        Volume old_vol(old_v);
        Volume new_vol(new_v);
        if ( sd.isValid() && old_vol.isSensitive() )   {
          new_vol.setSensitiveDetector(sd);
        }
        if ( c == geo_volume_t::Class() )  {
          Volume::Object *new_e, *old_e = (Volume::Object*)_data(old_vol);
          old_e->reflected = new_v;
          new_v->SetUserExtension(new_e = new Volume::Object(*old_e)); 
          new_e->reflected = old_v;
        }
        else if  ( c == geo_assembly_t::Class() )  {
          Assembly::Object *new_e, *old_e = (Assembly::Object*)_data(old_vol);
          old_e->reflected = new_v;
          new_v->SetUserExtension(new_e = new Assembly::Object(*old_e));
          new_e->reflected = old_v;
        }
        else if  ( c == TGeoVolumeMulti::Class() )   {
          VolumeMulti::Object *new_e, *old_e = (VolumeMulti::Object*)_data(old_vol);
          TGeoVolumeMulti* new_mv = (TGeoVolumeMulti*)new_v;
          TGeoVolumeMulti* old_mv = (TGeoVolumeMulti*)old_v;
          new_v->SetUserExtension(new_e = new VolumeMulti::Object(*old_e));
          old_e->reflected = new_v;
          new_e->reflected = old_v;
          for(int i=0, n=new_mv->GetNvolumes(); i<n; ++i)  {
            TGeoVolume* vol = new_mv->GetVolume(i);
            this->setShapeTitle(vol);
            (*this)(vol, old_mv->GetVolume(i), sd, set_bit);
          }
        }
        else
          except("dd4hep","VolumeImport: Unknown TGeoVolume sub-class: %s",new_v->IsA()->GetName());
 
        new_vol.setSensitiveDetector(old_vol.sensitiveDetector());
        new_vol.setVisAttributes(old_vol.visAttributes());
        new_vol.setLimitSet(old_vol.limitSet());
        new_vol.setRegion(old_vol.region());
 
        if ( set_bit >= 0 ) new_vol.setFlagBit(set_bit);
        for(Int_t i=0; i<new_v->GetNdaughters(); ++i)  {
          geo_node_t* pv = new_v->GetNode(i);
          geo_node_t* ov = old_v->GetNode(i);
          if ( !pv->GetUserExtension() )   {
            auto* e = (PlacedVolume::Object*)ov->geo_node_t::GetUserExtension();
            pv->geo_node_t::SetUserExtension(new PlacedVolume::Object(*e));
          }
          (*this)(pv->GetVolume(), ov->GetVolume(), sd, set_bit);
        }
      }
    }
  };
 
  TGeoVolume* MakeReflection(TGeoVolume* v, const char* newname=0)  {
    static TMap map(100);
    TGeoVolume* vol = (TGeoVolume*)map.GetValue(v);
    if ( vol ) {
      if (newname && newname[0]) v->SetName(newname);
      return vol;
    }
    vol = v->CloneVolume();
    if (!vol) {
      printout(ERROR,"MakeReflection", "Cannot clone volume %s\n", v->GetName());
      return nullptr;
    }
    map.Add(v, vol);
    std::string nam;
    if (newname && newname[0])  {
      nam = newname;
      vol->SetName(newname);
    }
    else   {
      nam = v->GetName();
      vol->SetName((nam+"_refl").c_str());
    }
    delete vol->GetNodes();
    vol->SetNodes(nullptr);
    vol->SetBit(TGeoVolume::kVolumeImportNodes, kFALSE);
    v->CloneNodesAndConnect(vol);
    // The volume is now properly cloned, but with the same shape.
    // Reflect the shape (if any) and connect it.
    if (v->GetShape())   {
      TGeoScale* scale = new TGeoScale( 1., 1.,-1.);
      TGeoShape* reflected_shape =
        TGeoScaledShape::MakeScaledShape((nam+"_shape_refl").c_str(), v->GetShape(), scale);
      vol->SetShape(reflected_shape);
      //vol->SetShape(v->GetShape());
    }
    // Reflect the daughters.
    Int_t nd = vol->GetNdaughters();
    if ( !nd ) return vol;
    TGeoVolume *new_vol;
    if ( !vol->GetFinder() ) {
      for (Int_t i=0; i < nd; i++) {
        TGeoNodeMatrix* node  = (TGeoNodeMatrix*)vol->GetNode(i);
        TGeoMatrix*     local = node->GetMatrix();
        //         printf("%s before\n", node->GetName());
        //         local->Print();
        Bool_t      reflected = local->IsReflection();
        TGeoMatrix* local_cloned = new TGeoCombiTrans(*local);
        local_cloned->RegisterYourself();
        node->SetMatrix(local_cloned);
        if (!reflected) {
          // We need to reflect only the translation and propagate to daughters.
          // H' = Sz * H * Sz
          local_cloned->ReflectZ(kTRUE);
          local_cloned->ReflectZ(kFALSE);
          //            printf("%s after\n", node->GetName());
          //            node->GetMatrix()->Print();
          new_vol = MakeReflection(node->GetVolume());
          node->SetVolume(new_vol);
          continue;
        }
        // The next daughter is already reflected, so reflect on Z everything and stop
        local_cloned->ReflectZ(kTRUE); // rot + tr
        //         printf("%s already reflected... After:\n", node->GetName());
        //         node->GetMatrix()->Print();
      }
      if ( vol->GetVoxels() ) vol->GetVoxels()->Voxelize();
      return vol;
    }
    // Volume is divided, so we have to reflect the division.
    //   printf("   ... divided %s\n", fFinder->ClassName());
    TGeoPatternFinder *new_finder = v->GetFinder()->MakeCopy(kTRUE);
    if (!new_finder) {
      printout(ERROR,"MakeReflection", "Could not copy finder for volume %s", v->GetName());
      return nullptr;
    }
    new_finder->SetVolume(vol);
    vol->SetFinder(new_finder);
    TGeoNodeOffset *nodeoff;
    new_vol = 0;
    for (Int_t i=0; i<nd; i++) {
      nodeoff = (TGeoNodeOffset*)vol->GetNode(i);
      nodeoff->SetFinder(new_finder);
      new_vol = MakeReflection(nodeoff->GetVolume());
      nodeoff->SetVolume(new_vol);
    }
    return vol;
  }
 
  int get_copy_number(TGeoVolume* par)    {
    TObjArray* a = par ? par->GetNodes() : 0;
    int copy_nr = (a ? a->GetEntries() : 0);
    return copy_nr;
  }
 
}
 
void ReflectionBuilder::execute()  const   {
  TGeoIterator next(detector.manager().GetTopVolume());
  bool print_active = isActivePrintLevel(DEBUG);
  TGeoNode *node;
  while ( (node=next()) ) {
    TGeoMatrix* matrix = node->GetMatrix();
    if (matrix->IsReflection()) {
      if ( print_active )  {
        printout(INFO,"ReflectionBuilder","Reflection matrix:");
        matrix->Print();
      }
      Volume vol(node->GetVolume());
      TGeoMatrix* mclone = new TGeoCombiTrans(*matrix);
      mclone->RegisterYourself();
      // Reflect just the rotation component
      //mclone->ReflectZ(kFALSE, kTRUE);
      if ( print_active )  {
        printout(INFO,"ReflectionBuilder","CLONE matrix:");
        mclone->Print();
      }
      TGeoNodeMatrix* nodematrix = (TGeoNodeMatrix*)node;
      nodematrix->SetMatrix(mclone);
      if ( print_active )  {
        printout(INFO,"ReflectionBuilder","Reflected volume: %s ", vol.name());
      }
      Volume refl = vol.reflect(vol.sensitiveDetector());
      node->SetVolume(refl.ptr());
    }
  }
}
 
PlacedVolume::Processor::Processor()   {
}
 
PlacedVolume::Processor::~Processor()   {
}
 
PlacedVolumeExtension::Parameterisation* PlacedVolumeExtension::Parameterisation::addref()  {
  ++refCount;
  return this;
}
void PlacedVolumeExtension::Parameterisation::release()  {
  --refCount;
  if ( 0 == refCount ) delete this;
}
 
PlacedVolumeExtension::PlacedVolumeExtension()
  : TGeoExtension(), volIDs()
{
  magic = magic_word();
  INCREMENT_COUNTER;
}
 
PlacedVolumeExtension::PlacedVolumeExtension(PlacedVolumeExtension&& c)
  : TGeoExtension(c), magic(std::move(c.magic)), refCount(0), volIDs(std::move(c.volIDs)) {
  INCREMENT_COUNTER;
}
 
PlacedVolumeExtension::PlacedVolumeExtension(const PlacedVolumeExtension& c)
  : TGeoExtension(), magic(c.magic), refCount(0), volIDs() {
  INCREMENT_COUNTER;
  volIDs = c.volIDs;
}
 
PlacedVolumeExtension::~PlacedVolumeExtension() {
  if ( this->params ) this->params->release();
  DECREMENT_COUNTER;
}
 
PlacedVolumeExtension& PlacedVolumeExtension::operator=(PlacedVolumeExtension&& copy)  {
  magic  = std::move(copy.magic);
  params = std::move(copy.params);
  volIDs = std::move(copy.volIDs);
  return *this;
}
PlacedVolumeExtension& PlacedVolumeExtension::operator=(const PlacedVolumeExtension& copy) {
  magic  = copy.magic;
  params = copy.params;
  volIDs = copy.volIDs;
  return *this;
}
 
TGeoExtension* PlacedVolumeExtension::Grab()   {
  ++this->refCount;
#ifdef ___print_vols
  else  cout << "Placement grabbed....." << endl;
#endif
  return this;
}
 
void PlacedVolumeExtension::Release() const  {
  PlacedVolumeExtension* ext = const_cast<PlacedVolumeExtension*>(this);
  --ext->refCount;
  if ( 0 == ext->refCount ) delete ext;
}
 
std::vector<PlacedVolumeExtension::VolID>::const_iterator
PlacedVolumeExtension::VolIDs::find(const std::string& name) const {
  for (Base::const_iterator i = this->Base::begin(); i != this->Base::end(); ++i)
    if (name == (*i).first)
      return i;
  return this->end();
}
 
std::pair<std::vector<PlacedVolumeExtension::VolID>::iterator, bool>
PlacedVolumeExtension::VolIDs::insert(const std::string& name, int value) {
  Base::iterator i = this->Base::begin();
  for (; i != this->Base::end(); ++i)
    if (name == (*i).first)
      break;
  //
  if (i != this->Base::end()) {
    return make_pair(i, false);
  }
  i = this->Base::emplace(this->Base::end(), name, value);
  return make_pair(i, true);
}
 
std::string PlacedVolumeExtension::VolIDs::str()  const   {
  std::stringstream str;
  str << std::hex;
  for(const auto& i : *this )   {
    str << i.first << "=" << std::setw(4) << std::right
        << std::setfill('0') << i.second << std::setfill(' ') << " ";
  }
  return str.str();
}
 
PlacedVolume::Object* PlacedVolume::data() const   {
  PlacedVolume::Object* o = _userExtension(*this);
  return o;
}
 
const char* PlacedVolume::type() const   {
  return m_element ? m_element->IsA()->GetName() : "UNKNOWN-PlacedVolume";
}
 
int PlacedVolume::copyNumber() const   {
  return m_element ? m_element->GetNumber() : -1;
}
 
Material PlacedVolume::material() const {
  return Material(m_element ? m_element->GetMedium() : 0);
}
 
Volume PlacedVolume::volume() const {
  return Volume(m_element ? m_element->GetVolume() : 0);
}
 
Volume PlacedVolume::motherVol() const {
  return Volume(m_element ? m_element->GetMotherVolume() : 0);
}
 
std::size_t PlacedVolume::num_daughters()  const   {
  return m_element ? m_element->GetNdaughters() : 0;
}
 
PlacedVolume PlacedVolume::daughter(std::size_t which)  const   {
  if ( m_element )  {
    if ( which < (std::size_t)m_element->GetNdaughters() )   {
      return m_element->GetDaughter(which);
    }
    except("Volume","+++ Access daughter %ld of %s [Has only %d daughters]",
       which, m_element->GetName(), m_element->GetNdaughters());
  }
  except("Volume","+++ Cannot access daughters of a non-existing volume!");
  return nullptr;
}
 
const PlacedVolume::VolIDs& PlacedVolume::volIDs() const {
  return _data(*this)->volIDs;
}
 
PlacedVolume& PlacedVolume::addPhysVolID(const std::string& nam, int value) {
  auto* o = _data(*this);
  if ( !o->params )   {
    o->volIDs.emplace_back(nam, value);
    return *this;
  }
  if ( value > 0 )    {
    except("PlacedVolume",
       "+++ addPhysVolID(%s): parameterised volumes only accept '0' is volID."
       "These automatically get overwritten with the copy number!",
       ptr()->GetName());
  }
  if ( !o->volIDs.empty() )    {
    except("PlacedVolume",
       "+++ addPhysVolID(%s): parameterised volumes can only host 1 physical volume ID."
       " vol id '%s' is already defined!", ptr()->GetName(), o->volIDs[0].first.c_str());
  }
  for(PlacedVolume pv : o->params->placements)  {
    auto* p = _data(pv);
    p->volIDs.emplace_back(nam, pv->GetNumber());
  }
  return *this;
}
 
const TGeoMatrix& PlacedVolume::matrix()  const    {
  if ( !isValid() )  {
    except("PlacedVolume","+++ matrix: Failed to access invalid PlacedVolume! [Invalid handle]");
  }
  return *(m_element->GetMatrix());  
}
 
Position PlacedVolume::position()  const    {
  const double* ptr = matrix().GetTranslation();
  return Position(ptr[0],ptr[1],ptr[2]);
}
 
std::string PlacedVolume::toString() const {
  std::stringstream str;
  Object* obj = _data(*this);
  str << m_element->GetName() << ":  vol='" << m_element->GetVolume()->GetName()
      << "' mat:'" << m_element->GetMatrix()->GetName()
      << "' volID[" << obj->volIDs.size() << "] ";
  for (VolIDs::const_iterator i = obj->volIDs.begin(); i != obj->volIDs.end(); ++i)
    str << (*i).first << "=" << (*i).second << "  ";
  str << std::ends;
  return str.str();
}
 
VolumeExtension::VolumeExtension()
: TGeoExtension(), region(), limits(), vis(), sens_det()   {
  INCREMENT_COUNTER;
}
 
VolumeExtension::~VolumeExtension() {
  detail::deletePtr(properties);
  region.clear();
  limits.clear();
  vis.clear();
  sens_det.clear();
  DECREMENT_COUNTER;
}
 
VolumeExtension::VolumeExtension(const VolumeExtension& c)   {
  this->copy(c);
}
 
VolumeExtension& VolumeExtension::operator=(const VolumeExtension& c)   {
  this->copy(c);
  return *this;
}
 
void VolumeExtension::copy(const VolumeExtension& c) {
  if ( this != &c )   {
    magic      = c.magic;
    region     = c.region;
    limits     = c.limits;
    vis        = c.vis;
    sens_det   = c.sens_det;
    referenced = c.referenced;
    detail::deletePtr(properties);
    if ( c.properties )    {
      properties = new TList();
      properties->SetOwner();
      TIter next(properties);
      TNamed *property;
      while ((property = (TNamed*)next())) properties->Add(new TNamed(*property));
    }
  }
}
 
TGeoExtension* VolumeExtension::Grab()  {
  VolumeExtension* ext = const_cast<VolumeExtension*>(this);
  ++ext->refCount;
#ifdef ___print_vols
  if ( ext->sens_det.isValid() )
    cout << "Volume grabbed with valid sensitive detector....." << endl;
  else
    cout << "Volume grabbed....." << endl;
#endif
  return ext;
}
 
void VolumeExtension::Release() const  {
  VolumeExtension* ext = const_cast<VolumeExtension*>(this);
  --ext->refCount;
  if ( 0 == ext->refCount )   {
#ifdef ___print_vols
    cout << "Volume deleted....." << endl;
#endif
    delete ext;
  }
  else  {
#ifdef ___print_vols
    cout << "VolumeExtension::Release::refCount:" << ext->refCount << endl;
#endif
  }
}
 
Volume::Volume(const std::string& nam) {
  m_element = _createTGeoVolume(nam,0,0);
}
 
Volume::Volume(const std::string& nam, const std::string& title) {
  m_element = _createTGeoVolume(nam,0,0);
  m_element->SetTitle(title.c_str());
}
 
Volume::Volume(const std::string& nam, const Solid& sol, const Material& mat) {
  m_element = _createTGeoVolume(nam, sol.ptr(), mat.ptr());
}
 
Volume::Volume(const std::string& nam, const std::string& title, const Solid& sol, const Material& mat) {
  m_element = _createTGeoVolume(nam, sol.ptr(), mat.ptr());
  m_element->SetTitle(title.c_str());
}
 
void Volume::enableCopyNumberCheck(bool value)    {
  s_verifyCopyNumbers = value;
}
 
Volume::Object* Volume::data() const   {
  Volume::Object* o = _userExtension(*this);
  return o;
}
 
const char* Volume::type() const   {
  return m_element ? m_element->IsA()->GetName() : "UNKNOWN-Volume";
}
 
Volume Volume::reflect()  const   {
  return this->reflect(this->sensitiveDetector());
}
    
Volume Volume::reflect(SensitiveDetector sd)  const   {
  if ( m_element )   {
    VolumeImport imp;
    Object* o = data();
    if ( !o->reflected.isValid() )  {
      TGeoVolume* vol = MakeReflection(m_element);
      imp(vol, m_element, sd, Volume::REFLECTED);
      o->reflected = vol;
    }
    return o->reflected;
  }
  except("dd4hep","Volume: Attempt to reflect an invalid Volume handle.");
  return *this;
}
 
Volume& Volume::import()    {
  if ( m_element )   {
    VolumeImport imp;
    imp(m_element);
    return *this;
  }
  except("dd4hep","Volume: Attempt to import an invalid Volume handle.");
  return *this;
}
 
void Volume::setFlagBit(unsigned int bit)   {
  if ( bit <= 31 )   {
    data()->flags |= 1<<bit;
    return;
  }
  except("Volume","+++ Volume flag bit outsize range [0...31]: %d",bit);
}
 
bool Volume::testFlagBit(unsigned int bit)   const    {
  if ( bit <= 31 )   {
    return (data()->flags & 1<<bit) != 0;
  }
  except("Volume","+++ Volume flag bit outsize range [0...31]: %d",bit);
  return false; // Anyhow never called. Only to satisfy the compiler.
}    
 
bool Volume::isReflected()   const    {
  return testFlagBit(REFLECTED);
}
 
bool Volume::isAssembly()   const   {
  return m_element ? m_element->IsAssembly() : false;
}    
 
unsigned char Volume::setSmartlessValue(unsigned char new_value)  {
  Object* obj = _data(*this);
  unsigned char tmp = obj->smartLess;
  obj->smartLess = new_value;
  return tmp;
}
 
unsigned char Volume::smartlessValue()  const  {
  return _data(*this)->smartLess;
}
 
Volume Volume::divide(const std::string& divname, int iaxis, int ndiv,
                      double start, double step, int numed, const char* option)   {
  TGeoVolume* p = m_element;
  if ( p )  {
    TGeoVolume* mvp = p->Divide(divname.c_str(), iaxis, ndiv, start, step, numed, option);
    if ( mvp )   {
      VolumeImport imp;
      imp(mvp);
      return VolumeMulti(mvp);
    }
    except("dd4hep","Volume: Failed to divide volume %s -> %s [Invalid result]",
           p->GetName(), divname.c_str());
  }
  except("dd4hep","Volume: Attempt to divide an invalid logical volume.");
  return nullptr;
}
 
PlacedVolume _addNode(TGeoVolume* par, TGeoVolume* daughter, int id, TGeoMatrix* transform) {
  TGeoVolume* parent = par;
  if ( !parent )   {
    except("dd4hep","Volume: Attempt to assign daughters to an invalid physical parent volume.");
  }
  else if ( !daughter )   {
    except("dd4hep","Volume: Attempt to assign an invalid physical daughter volume.");
  }
  else if ( !transform )   {
    except("dd4hep","Volume: Attempt to place volume without placement matrix.");
  }
  if ( transform != detail::matrix::_identity() ) {
    std::string nam = std::string(daughter->GetName()) + "_placement";
    transform->SetName(nam.c_str());
  }
  TGeoShape* shape = daughter->GetShape();
  // Need to fix the daughter's BBox of assemblies, if the BBox was not calculated....
  if ( shape->IsA() == TGeoShapeAssembly::Class() )  {
    TGeoShapeAssembly* as = (TGeoShapeAssembly*)shape;
    as->NeedsBBoxRecompute();
    as->ComputeBBox();
  }
  const Double_t* r = transform->GetRotationMatrix();
  if ( r )   {
    Double_t test_rot = r[0] + r[4] + r[8] - 3.0;
    if ( TMath::Abs(test_rot) < s_rotation_test_limit )
      transform->ResetBit(TGeoMatrix::kGeoRotation);
    else
      transform->SetBit(TGeoMatrix::kGeoRotation);
 
    if ( transform->IsRotation() )   {
      Double_t det =
        r[0]*r[4]*r[8] + r[3]*r[7]*r[2] + r[6]*r[1]*r[5] -
        r[2]*r[4]*r[6] - r[5]*r[7]*r[0] - r[8]*r[1]*r[3];
      if ( det < 0e0 )   {
        transform->SetBit(TGeoMatrix::kGeoReflection);
        printout(DEBUG, "PlacedVolume",
                 "REFLECTION: (x.Cross(y)).Dot(z): %8.3g Parent: %s [%s] Daughter: %s [%s]",
                 det, par->GetName(), par->IsA()->GetName(),
                 daughter->GetName(), daughter->IsA()->GetName());
      }
    }
  }
  geo_node_t* n {nullptr};
  TString nam_id = TString::Format("%s_%d", daughter->GetName(), id);
  if ( s_verifyCopyNumbers )   {
    n = static_cast<geo_node_t*>(parent->GetNode(nam_id));
    if ( n != 0 )  {
      printout(ERROR,"PlacedVolume","++ Attempt to place already exiting node %s",(const char*)nam_id);
    }
  }
  /* n = */ parent->AddNode(daughter, id, transform);
  //n = static_cast<geo_node_t*>(parent->GetNode(nam_id));
  n = static_cast<geo_node_t*>(parent->GetNodes()->Last());
  if ( nam_id != n->GetName() )   {
    printout(ERROR,"PlacedVolume","++ FAILED to place node %s",(const char*)nam_id);
  }
  PlacedVolume::Object* extension = new PlacedVolume::Object();
  n->geo_node_t::SetUserExtension(extension);
  return PlacedVolume(n);
}
 
PlacedVolume _addNode(TGeoVolume* par, Volume daughter, int copy_nr, const Rotation3D& rot3D)   {
  TGeoRotation r;
  double elements[9];
  rot3D.GetComponents(elements);
  r.SetMatrix(elements);
  auto matrix = std::make_unique<TGeoCombiTrans>(TGeoTranslation(0,0,0),r);
  return _addNode(par, daughter, copy_nr, matrix.release());
}
 
PlacedVolume _addNode(TGeoVolume* par, Volume daughter, int copy_nr, const Transform3D& tr)   {
  TGeoRotation r;
  double elements[9];
  Position   pos3D;
  Rotation3D rot3D;
  tr.GetRotation(rot3D);
  tr.GetTranslation(pos3D);
  rot3D.GetComponents(elements);
  r.SetMatrix(elements);
  auto matrix = std::make_unique<TGeoCombiTrans>(TGeoTranslation(pos3D.x(), pos3D.y(), pos3D.z()),r);
  return _addNode(par, daughter, copy_nr, matrix.release());
}
 
PlacedVolume Volume::placeVolume(const Volume& volume, TGeoMatrix* tr) const    {
  return _addNode(m_element, volume, get_copy_number(m_element), tr);
}
 
PlacedVolume Volume::placeVolume(const Volume& volume, int copy_nr, TGeoMatrix* tr) const    {
  return _addNode(m_element, volume, copy_nr, tr);
}
 
PlacedVolume Volume::placeVolume(const Volume& volume, const Transform3D& trans) const {
  return _addNode(m_element, volume, get_copy_number(m_element), trans);
}
 
PlacedVolume Volume::placeVolume(const Volume& volume) const {
  return _addNode(m_element, volume, get_copy_number(m_element), detail::matrix::_identity());
}
 
PlacedVolume Volume::placeVolume(const Volume& volume, const Position& pos) const {
  return _addNode(m_element, volume, get_copy_number(m_element), detail::matrix::_translation(pos));
}
 
PlacedVolume Volume::placeVolume(const Volume& volume, const RotationZYX& rot) const {
  return _addNode(m_element, volume, get_copy_number(m_element), detail::matrix::_rotationZYX(rot));
}
 
PlacedVolume Volume::placeVolume(const Volume& volume, const Rotation3D& rot) const {
  return _addNode(m_element, volume, get_copy_number(m_element), rot);
}
 
PlacedVolume Volume::placeVolume(const Volume& volume, int copy_no, const Transform3D& trans) const {
  return _addNode(m_element, volume, copy_no, trans);
}
 
PlacedVolume Volume::placeVolume(const Volume& volume, int copy_no) const {
  return _addNode(m_element, volume, copy_no, detail::matrix::_identity());
}
 
PlacedVolume Volume::placeVolume(const Volume& volume, int copy_no, const Position& pos) const {
  return _addNode(m_element, volume, copy_no, detail::matrix::_translation(pos));
}
 
PlacedVolume Volume::placeVolume(const Volume& volume, int copy_no, const RotationZYX& rot) const {
  return _addNode(m_element, volume, copy_no, detail::matrix::_rotationZYX(rot));
}
 
PlacedVolume Volume::placeVolume(const Volume& volume, int copy_no, const Rotation3D& rot) const {
  return _addNode(m_element, volume, copy_no, rot);
}
 
PlacedVolume Volume::replicate(const Volume entity, ReplicationAxis axis,
                   size_t count, double inc, double start)
{
  Transform3D offset(1e0,  0e0,  0e0,  axis == X_axis ? start : 0e0,
             0e0,  1e0,  0e0,  axis == Y_axis ? start : 0e0,
             0e0,  0e0,  1e0,  axis == Z_axis ? start : 0e0);
  Transform3D tr(1e0,  0e0,  0e0,  axis == X_axis ? inc : 0e0,
         0e0,  1e0,  0e0,  axis == Y_axis ? inc : 0e0,
         0e0,  0e0,  1e0,  axis == Z_axis ? inc : 0e0);
  PlacedVolume pv = paramVolume1D(offset, entity, count, tr);
  auto* data = pv.data();
  data->params->flags = axis | REPLICATED;
  return pv;
}
 
PlacedVolume Volume::paramVolume1D(Volume entity, size_t count, const Position& inc)   {
  return paramVolume1D(Transform3D(), entity, count, Transform3D(inc));
}
 
PlacedVolume Volume::paramVolume1D(Volume entity, size_t count, const RotationZYX& inc)  {
  return paramVolume1D(Transform3D(), entity, count, Transform3D(inc));
}
 
PlacedVolume Volume::paramVolume1D(Volume entity, size_t count, const Transform3D& inc)   {
  return paramVolume1D(Transform3D(), entity, count, inc);
}
 
PlacedVolume Volume::paramVolume1D(const Transform3D& start,
                   Volume entity,
                   size_t count,
                   const Transform3D& trafo)
{
  Transform3D tr(start);
  PlacedVolume pv =
    _addNode(m_element, entity, get_copy_number(m_element), detail::matrix::_transform(tr));
  auto* data = pv.data();
  if ( pv->GetNdaughters() > 1 )   {
    except("Volume","paramVolume1D: Mother %s has too many daughters: %ld "
       "Parameterized volumes may only have one single daughter!",
       ptr()->GetName(), pv.volume()->GetNdaughters());
  }
  data->params = new PlacedVolumeExtension::Parameterisation();
  data->params->addref();
  data->params->flags = PARAMETERIZED;
  data->params->start = start;
  data->params->trafo1D.first  = trafo;
  data->params->trafo1D.second = count;
  data->params->trafo2D.second = 0;
  data->params->trafo3D.second = 0;
  data->params->placements.emplace_back(pv);
  for(size_t i=1; i < count; ++i)    {
    tr *= trafo;
    PlacedVolume ppv =
      _addNode(m_element, entity, get_copy_number(m_element), detail::matrix::_transform(tr));
    data->params->placements.emplace_back(ppv);
    ppv.data()->params = data->params->addref();
  }
  return pv;
}
 
PlacedVolume Volume::paramVolume2D(Volume entity,
                   size_t count_1,
                   const Transform3D& trafo_1,
                   size_t count_2,
                   const Transform3D& trafo_2)
{
  return paramVolume2D(Transform3D(), entity, count_1, trafo_1, count_2, trafo_2);
}
 
PlacedVolume Volume::paramVolume2D(const Transform3D& start,
                   Volume entity,
                   size_t count_1,
                   const Position& pos_1,
                   size_t count_2,
                   const Position& pos_2)
{
  return paramVolume2D(start, entity, count_1, Transform3D(pos_1), count_2, Transform3D(pos_2));
}
 
PlacedVolume Volume::paramVolume2D(Volume entity,
                   size_t count_1,
                   const Position& pos_1,
                   size_t count_2,
                   const Position& pos_2)
{
  return paramVolume2D(Transform3D(), entity, count_1, Transform3D(pos_1), count_2, Transform3D(pos_2));
}
 
PlacedVolume Volume::paramVolume2D(const Transform3D& start,
                   Volume entity,
                   size_t count_1,
                   const Transform3D& trafo_1,
                   size_t count_2,
                   const Transform3D& trafo_2)
{
  PlacedVolume pv =
    _addNode(m_element, entity, get_copy_number(m_element), detail::matrix::_transform(start));
  auto* data = pv.data();
  if ( pv->GetNdaughters() > 1 )   {
    except("Volume","paramVolume1D: Mother %s has too many daughters: %ld "
       "Parameterized volumes may only have one single daughter!",
       ptr()->GetName(), pv.volume()->GetNdaughters());
  }
  data->params = new PlacedVolumeExtension::Parameterisation();
  data->params->addref();
  data->params->flags = PARAMETERIZED;
  data->params->start = start;
  data->params->trafo1D.first  = trafo_1;
  data->params->trafo1D.second = count_1;
  data->params->trafo2D.first  = trafo_2;
  data->params->trafo2D.second = count_2;
  data->params->trafo3D.second = 0;
  data->params->placements.emplace_back(pv);
  Transform3D tr2(start);
  for(size_t j=0; j < count_2; ++j)    {
    Transform3D tr1 = tr2;
    for(size_t i = 0; i < count_1; ++i)    {
      if ( !( i == 0 && j == 0 ) )   {
    PlacedVolume ppv =
      _addNode(m_element, entity, get_copy_number(m_element), detail::matrix::_transform(tr1));
    data->params->placements.emplace_back(ppv);
    ppv.data()->params = data->params->addref();
      }
      tr1 *= trafo_1;
    }
    tr2 *= trafo_2;
  }
  return pv;
}
 
PlacedVolume Volume::paramVolume3D(const Transform3D& start,
                   Volume entity,
                   size_t count_1,
                   const Position& pos_1,
                   size_t count_2,
                   const Position& pos_2,
                   size_t count_3,
                   const Position& pos_3)
{
  return paramVolume3D(start, entity, 
               count_1, Transform3D(pos_1),
               count_2, Transform3D(pos_2),
               count_3, Transform3D(pos_3));
}
 
PlacedVolume Volume::paramVolume3D(Volume entity,
                   size_t count_1,
                   const Position& pos_1,
                   size_t count_2,
                   const Position& pos_2,
                   size_t count_3,
                   const Position& pos_3)
{
  return paramVolume3D(Transform3D(), entity,
               count_1, Transform3D(pos_1),
               count_2, Transform3D(pos_2),
               count_3, Transform3D(pos_3));
}
 
PlacedVolume Volume::paramVolume3D(const Transform3D& start,
                   Volume entity,
                   size_t count_1,
                   const Transform3D& trafo_1,
                   size_t count_2,
                   const Transform3D& trafo_2,
                   size_t count_3,
                   const Transform3D& trafo_3)
{
  PlacedVolume pv =
    _addNode(m_element, entity, get_copy_number(m_element), detail::matrix::_transform(start));
  auto* data = pv.data();
  if ( pv->GetNdaughters() > 1 )   {
    except("Volume","paramVolume1D: Mother %s has too many daughters: %ld "
       "Parameterized volumes may only have one single daughter!",
       ptr()->GetName(), pv.volume()->GetNdaughters());
  }
  data->params = new PlacedVolumeExtension::Parameterisation();
  data->params->addref();
  data->params->flags = PARAMETERIZED;
  data->params->start = start;
  data->params->trafo1D.first  = trafo_1;
  data->params->trafo1D.second = count_1;
  data->params->trafo2D.first  = trafo_2;
  data->params->trafo2D.second = count_2;
  data->params->trafo3D.first  = trafo_3;
  data->params->trafo3D.second = count_3;
  data->params->placements.emplace_back(pv);
  Transform3D tr3(start);
  for(size_t k=0; k < count_3; ++k)    {
    Transform3D tr2 = tr3;
    for(size_t j=0; j < count_2; ++j)    {
      Transform3D tr1 = tr2;
      for(size_t i = 0; i < count_1; ++i)    {
    if ( !( i == 0 && j == 0 && k == 0 ) )   {
      PlacedVolume ppv =
        _addNode(m_element, entity, get_copy_number(m_element), detail::matrix::_transform(tr1));
      data->params->placements.emplace_back(ppv);
      ppv.data()->params = data->params->addref();
    }
    tr1 *= trafo_1;
      }
      tr2 *= trafo_2;
    }
    tr3 *= trafo_3;
  }
  return pv;
}
 
const Volume& Volume::setOption(const std::string& opt) const {
  if ( isValid() )   {
    m_element->SetOption(opt.c_str());
    return *this;
  }
  throw std::runtime_error("dd4hep: Attempt to access invalid handle of type: PlacedVolume");
}
 
std::string Volume::option() const {
  return m_element->GetOption();
}
 
const Volume& Volume::setMaterial(const Material& mat) const {
  if (mat.isValid()) {
    TGeoMedium* medium = mat._ptr<TGeoMedium>();
    if (medium) {
      m_element->SetMedium(medium);
      return *this;
    }
    throw std::runtime_error("dd4hep: Volume: Medium " + std::string(mat.name()) + " is not registered with geometry manager.");
  }
  throw std::runtime_error("dd4hep: Volume: Attempt to assign invalid material.");
}
 
Material Volume::material() const {
  return Material(m_element->GetMedium());
}
 
const Volume& Volume::setVisAttributes(const VisAttr& attr) const {
  if ( attr.isValid() ) {
    VisAttr::Object* vis = attr.data<VisAttr::Object>();
    TColor* col = vis->color;
    if ( col )   {
      int draw_style = vis->drawingStyle;
      int line_style = vis->lineStyle;
      int col_num    = col->GetNumber();
      int col_tr_num = vis->colortr->GetNumber();
      m_element->SetVisibility(vis->visible ? kTRUE : kFALSE);
      m_element->SetVisContainers(kTRUE);
      m_element->SetVisDaughters(vis->showDaughters ? kTRUE : kFALSE);
      printout(DEBUG,"setVisAttributes",
               "Set color %3d transparent(alpha:%.3f): %3d [%02X,%02X,%02X] DrawingStyle:%9s LineStyle:%6s for volume %s",
               col_num, vis->alpha, col_tr_num,
               int(255*col->GetRed()),
               int(255*col->GetGreen()),
               int(255*col->GetBlue()),
               draw_style == VisAttr::SOLID ? "Solid" : "Wireframe",
               line_style == VisAttr::SOLID ? "Solid" : "Dashed",
               name()
               );
      m_element->SetLineWidth(10);
      m_element->SetLineColor(col_num);
      m_element->SetFillColor(col_tr_num);
      if (draw_style == VisAttr::SOLID) {
        m_element->SetFillStyle(1001);   // Root: solid
 
#if ROOT_VERSION_CODE >= ROOT_VERSION(6,29,0)
        // Set directly transparency to the volume, NOT to the material as for ROOT < 6.29
        m_element->ResetTransparency(Char_t((1.0-vis->alpha)*100));
#else
        // As suggested by Valentin Volkl https://sft.its.cern.ch/jira/browse/DDFORHEP-20
        //
        // According to https://root.cern.ch/phpBB3/viewtopic.php?t=2309#p66013
        // a transparency>50 will make a volume invisible in the normal pad.
        // Hence: possibly restrict transparency to a maximum of 50.
        //        but let's see first how this behaves.
        m_element->SetTransparency(Char_t((1.0-vis->alpha)*100));
#endif
      }
      else {
        printout(DEBUG,"setVisAttributes","Set to wireframe vis:%s",name());
        m_element->SetLineColor(kBlack);
        m_element->SetFillColor(0);
        m_element->SetFillStyle(0);      // Root: hollow
      }
      if (line_style == VisAttr::SOLID)  // Root line style: 1=solid, 2=dash, 3=dot, 4=dash-dot.
        m_element->SetLineStyle(1);
      else if (line_style == VisAttr::DASHED)
        m_element->SetLineStyle(2);
      else
        m_element->SetLineStyle(line_style);
    }
    else   {
      except("Volume","setVisAttributes: encountered valid, but badly initialized visattr: %s",attr.name());
    }
  }
  Volume::Object* o = _userExtension(*this);
  if ( o ) o->vis = attr;
  return *this;
}
 
const Volume& Volume::setVisAttributes(const Detector& description, const std::string& nam) const {
  if (!nam.empty()) {
    VisAttr attr = description.visAttributes(nam);
    setVisAttributes(attr);
  }
  return *this;
}
 
const Volume& Volume::setAttributes(const Detector& description, const std::string& rg, const std::string& ls, const std::string& vis) const {
  if (!rg.empty())
    setRegion(description.region(rg));
  if (!ls.empty())
    setLimitSet(description.limitSet(ls));
  setVisAttributes(description, vis);
  return *this;
}
 
VisAttr Volume::visAttributes() const {
  Object* o = _data(*this, false);
  if (o) return o->vis;
  return VisAttr();
}
 
const Volume& Volume::setSolid(const Solid& sol) const {
  m_element->SetShape(sol);
  return *this;
}
 
Solid Volume::solid() const {
  return Solid((*this)->GetShape());
}
 
Box Volume::boundingBox() const {
  Box box = this->solid();
  if ( box.isValid() )   {
    return box;
  }
  else if ( !isValid() )   {
    except("dd4hep","Volume: Cannot access the bounding box of an invalid volume [Invalid Handle]!");
  }
  except("dd4hep","Volume: Cannot access the bounding box an object of type: %s shape: %s",
         this->ptr()->IsA()->GetName(), this->ptr()->GetShape()->IsA()->GetName());
  return box;
}
 
const Volume& Volume::setRegion(const Detector& description, const std::string& nam) const {
  if (!nam.empty()) {
    return setRegion(description.region(nam));
  }
  return *this;
}
 
const Volume& Volume::setRegion(const Region& obj) const {
  _data(*this)->region = obj;
  return *this;
}
 
Region Volume::region() const {
  return _data(*this)->region;
}
 
const Volume& Volume::setLimitSet(const Detector& description, const std::string& nam) const {
  if (!nam.empty()) {
    return setLimitSet(description.limitSet(nam));
  }
  return *this;
}
 
const Volume& Volume::setLimitSet(const LimitSet& obj) const {
  _data(*this)->limits = obj;
  return *this;
}
 
LimitSet Volume::limitSet() const {
  return _data(*this)->limits;
}
 
const Volume& Volume::setSensitiveDetector(const SensitiveDetector& obj) const {
  //cout << "Setting sensitive detector '" << obj.name() << "' to volume:" << ptr() << " " << name() << endl;
  _data(*this)->sens_det = obj;
  return *this;
}
 
Ref_t Volume::sensitiveDetector() const {
  const Object* o = _data(*this);
  return o->sens_det;
}
 
bool Volume::isSensitive() const {
  return _data(*this)->sens_det.isValid();
}
 
bool Volume::hasProperties()  const   {
  return _data(*this)->properties != nullptr;
}
 
void Volume::addProperty(const std::string& nam, const std::string& val) const  {
  auto* o = _data(*this);
  if ( !o->properties )   {
    o->properties = new TList();
    o->properties->SetOwner();
  }
  TNamed *prop = (TNamed*)o->properties->FindObject(nam.c_str());
  if ( !prop )   {
    o->properties->Add(new TNamed(nam.c_str(), val.c_str()));
    return;
  }
  except("Volume::addProperty", "Volume: '%s' Property '%s' is already set!",
     ptr()->GetName(), nam.c_str());
}
 
std::string Volume::getProperty(const std::string& nam, const std::string& default_val)   const {
  const auto* o = _data(*this);
  if ( !o->properties )   {
    return default_val;
  }
  TNamed *prop = (TNamed*)o->properties->FindObject(nam.c_str());
  if ( prop ) return prop->GetTitle();
  return default_val;
}
 
Assembly::Assembly(const std::string& nam) {
  m_element = _createTGeoVolumeAssembly(nam);
}
 
VolumeMulti::VolumeMulti(const std::string& nam, Material mat) {
  m_element = _createTGeoVolumeMulti(nam, mat.ptr());
}
 
void VolumeMulti::verifyVolumeMulti()   {
  if ( m_element )  {
    // This will lead to an exception if the type is not TGeoVolumeMulti
    TGeoVolumeMulti* multi = detail::safe_cast<TGeoVolumeMulti>::cast(m_element);
    if ( multi )  {
      import();
      return;
    }
    // Force a bad cast exception
    Handle<TGeoVolumeMulti> handle(m_element);
    if ( handle.isValid() )  {}
  }
}
 
std::string dd4hep::toStringMesh(PlacedVolume place, int prec)   {
  Volume       vol   = place->GetVolume();
  TGeoMatrix*  mat   = place->GetMatrix();
  Solid        sol   = vol.solid();
  std::stringstream os;
  struct _numbering {
    double adjust(double value)  const   {
      if ( std::abs(value) < TGeoShape::Tolerance() )
        return 0.0;
      return value/dd4hep::cm;
    }
  } _vertex;
 
  if ( vol->IsA() == TGeoVolumeAssembly::Class() )    {
    for(int i=0; i<vol->GetNdaughters(); ++i)  {
      os << toStringMesh(vol->GetNode(i), prec) << std::endl;
    }
    return os.str();
  }
 
  // Prints shape parameters
  int nvert = 0, nsegs = 0, npols = 0;
  sol->GetMeshNumbers(nvert, nsegs, npols);
  Double_t* points = new Double_t[3*nvert];
  sol->SetPoints(points);
 
  os << std::setw(16) << std::left << sol->IsA()->GetName()
     << " " << nvert << " Mesh-points:" << std::endl;
  os << std::setw(16) << std::left << sol->IsA()->GetName() << " " << sol->GetName()
     << " N(mesh)=" << sol->GetNmeshVertices()
     << "  N(vert)=" << nvert << "  N(seg)=" << nsegs << "  N(pols)=" << npols << std::endl;
    
  for(int i=0; i<nvert; ++i)   {
    Double_t* p = points + 3*i;
    Double_t global[3], local[3] = {p[0], p[1], p[2]};
    mat->LocalToMaster(local, global);
    os << std::setw(16) << std::left << sol->IsA()->GetName() << " " << std::setw(3) << std::left << i
       << " Local  ("  << std::setw(7) << std::setprecision(prec) << std::fixed << std::right << _vertex.adjust(local[0])
       << ", "         << std::setw(7) << std::setprecision(prec) << std::fixed << std::right << _vertex.adjust(local[1])
       << ", "         << std::setw(7) << std::setprecision(prec) << std::fixed << std::right << _vertex.adjust(local[2])
       << ") Global (" << std::setw(7) << std::setprecision(prec) << std::fixed << std::right << _vertex.adjust(global[0])
       << ", "         << std::setw(7) << std::setprecision(prec) << std::fixed << std::right << _vertex.adjust(global[1])
       << ", "         << std::setw(7) << std::setprecision(prec) << std::fixed << std::right << _vertex.adjust(global[2])
       << ")" << std::endl;
  }
  Box box = sol;
  const Double_t* org = box->GetOrigin();
  os << std::setw(16) << std::left << sol->IsA()->GetName()
     << " Bounding box: "
     << " dx="        << std::setw(7) << std::setprecision(prec) << std::fixed << std::right << _vertex.adjust(box->GetDX())
     << " dy="        << std::setw(7) << std::setprecision(prec) << std::fixed << std::right << _vertex.adjust(box->GetDY())
     << " dz="        << std::setw(7) << std::setprecision(prec) << std::fixed << std::right << _vertex.adjust(box->GetDZ())
     << " Origin: x=" << std::setw(7) << std::setprecision(prec) << std::fixed << std::right << _vertex.adjust(org[0])
     << " y="         << std::setw(7) << std::setprecision(prec) << std::fixed << std::right << _vertex.adjust(org[1])
     << " z="         << std::setw(7) << std::setprecision(prec) << std::fixed << std::right << _vertex.adjust(org[2])
     << std::endl;
  
  delete [] points;
  return os.str();
}