GeoGeometry.cxx

Go to the documentation of this file.

//
// GeoGeometry
//
// GeoGeometry is single geometry constructed using the TGeoManager
//
// Author:  S. Kasahara 11/03
//
// Based on R. Hatcher's UgliGeometry/UgliGeometry class.

#include <iostream>
#include <cassert>
#include <math.h>
#include <set>
using namespace std;

#include <TGeoManager.h>
#include <TGeoMatrix.h>
#include <TGeoNode.h>
#include <TGeoBBox.h>
#include <TGeoArb8.h>
#include <TList.h>

#include "Conventions/Munits.h"
#include "Conventions/MinosMaterial.h"
#include "Util/UtilString.h"
#include "GeoGeometry/Geo.h"
#include "GeoGeometry/GeoGeometry.h"
#include "GeoGeometry/GeoMediumMap.h"
#include "GeoGeometry/GeoScintPlnVolume.h"
#include "GeoGeometry/GeoSteelPlnVolume.h"
#include "GeoGeometry/GeoStripVolume.h"
#include "GeoGeometry/GeoScintMdlVolume.h"
#include "GeoGeometry/GeoStripNode.h"
#include "GeoGeometry/GeoScintMdlNode.h"
#include "GeoGeometry/GeoScintPlnNode.h"
#include "GeoGeometry/GeoSteelPlnNode.h"
#include "GeoGeometry/GeoMedium.h"
#include "UgliGeometry/UgliDbiTables.h"
#include "UgliGeometry/UgliDbiStructHash.h"
#include "MessageService/MsgService.h"
#include "Plex/PlexVetoShieldHack.h"

ClassImp(GeoGeometry)

CVSID("$Id: GeoGeometry.cxx,v 1.87 2008/01/09 06:38:11 schubert Exp $");

//_____________________________________________________________________________
GeoGeometry::GeoGeometry() : fVldRange(),fAppType(Geo::kRecons),
                fScale(1),fHallPath(""),fGeoManager(0),fGeoShield(),
                fIsSwimMedia(false),fMediumMap(0),fRot45(0),fRot315(0) {
  // Default constructor, used for i/o

  MSG("Geo",Msg::kDebug) << "GeoGeometry def ctor @ " << this << endl;

  UpdateGlobalManager(); // this is okay, because fGeoManager ptr stored

  // fScale has been set before this step
  for ( int ism = 0; ism < 2; ism++ ) { 
     fSMZMin[ism] = +999.*fScale; 
     fSMZMax[ism] = -999.*fScale; 
  }

}

//_____________________________________________________________________________
GeoGeometry::GeoGeometry(const VldContext &vldc, Geo::EAppType apptype) : 
  fVldRange(),fAppType(apptype),fScale(1),fHallPath(""),
  fGeoManager(0),fGeoShield(this),fIsSwimMedia(false),fMediumMap(0),
  fRot45(0),fRot315(0) {
  // Normal constructor
  // Geometry is built for application type Geo::kRecons by default.
  // The application type determines the units the geometry will be built
  // in.  The default apptype = Geo::kRecons will build the geometry in
  // Minos standard units of meters.  Specifying apptype = Geo::kVMC will
  // build the geometry in TVirtualMC standard units of cm.
  
  MSG("Geo",Msg::kDebug) << "GeoGeometry normal ctor @ " << this << endl;

  UpdateGlobalManager(); // set gGeoManager to null before building new 

  MSG("Geo",Msg::kInfo) << "GeoGeometry build for validity " 
                        << vldc << "." << endl;
  
  // Silence TGeoManager info messages
  Int_t saveErrorIgnoreLevel = gErrorIgnoreLevel;
  gErrorIgnoreLevel = kWarning;
  
  //MsgStream* msgeo = MsgService::Instance() -> GetStream("Geo");
  //msgeo -> SetLogLevel(Msg::kVerbose);
  
  fScale = Geo::GetScale(apptype);

  // fScale must be set before this step
  for ( int ism = 0; ism < 2; ism++ ) { 
     fSMZMin[ism] = +999.*fScale; 
     fSMZMax[ism] = -999.*fScale; 
  }

  fGeoManager = new TGeoManager(vldc.AsString("s1ac-"),"a MINOS geometry");
  UpdateGlobalManager(); // set gGeoManager to fGeoManager

  fMediumMap = new GeoMediumMap(this,vldc);
  
  // Build all materials and detector geometry
  this -> BuildAll(vldc);
  
  fGeoManager->CloseGeometry();

  TGeoVolume* volume = fGeoManager->GetTopVolume();
  UpdateNodeMatrices(volume);
  // Invoke CdTop to position back to top (MARS) volume.
  fGeoManager->CdTop();
  
  //this -> DumpVolume("MARS");
  
  fGeoManager->SetVisLevel(4);
  fGeoManager->SetVisOption(0);

  MSG("Geo",Msg::kInfo) << "GeoGeometry build complete." << endl;
  gErrorIgnoreLevel = saveErrorIgnoreLevel;

}

//_____________________________________________________________________________
GeoGeometry::~GeoGeometry() {
  // delete all the owned sub-objects

  MSG("Geo",Msg::kDebug) << "GeoGeometry dtor @ " << this << endl;

  UpdateGlobalManager();

  fPlaneMap.clear();
  if ( fMediumMap ) delete fMediumMap; fMediumMap = 0;
  
  if ( fGeoManager ) {
    delete fGeoManager; fGeoManager = 0;
  }
 
  UpdateGlobalManager(); // set gGeoManager to null to avoid further use
  
}

//_____________________________________________________________________________
void GeoGeometry::Draw(Option_t* volname) {
  // Public. Draw this geometry from volume volname (default "HALL")

  UpdateGlobalManager();

  std::string volnamestr(volname);
  if (volnamestr.empty()) volnamestr = "HALL"; 
   
  MSG("Geo",Msg::kInfo) << "GeoGeometry::Draw() from top volume " 
                        << volnamestr.empty() << "." << endl;

  if ( !fGeoManager ) {
    MSG("Geo",Msg::kWarning) << "Null fGeoManager" << endl;
    return;
  }
  
  TGeoVolume* topVol = fGeoManager -> GetVolume(volnamestr.c_str());
  if ( !topVol ) {
    MSG("Geo",Msg::kWarning) << "No volume of name " << volnamestr.c_str() 
                             << endl;
    return;
  }
  //fGeoManager->SetTopVisible();
  topVol -> SetVisContainers();
  topVol -> Draw();
   
}

//_____________________________________________________________________________
void GeoGeometry::DumpVolume(std::string volname, std::string preface) const {
  // Print information about volume down through all daughter nodes

  UpdateGlobalManager();

  TGeoVolume* volume = fGeoManager->GetVolume(volname.c_str());
  if ( !volume ) {
    cout << "No volume of name \"" << volname.c_str() << "\"!" << endl;
    return;
  }
  
  cout << preface.c_str() << "Volume " << volume->GetName() << ", shape:" 
       << endl;
  cout << preface.c_str();
  volume->GetShape()->InspectShape();
  
  Int_t ndaughter = volume -> GetNdaughters();
  if ( ndaughter > 0 ) 
    cout << preface.c_str() << ndaughter << " daughter nodes:" << endl;
  else
    cout << preface.c_str() << "No daughter nodes." << endl;
  
  for ( int id = 0; id < ndaughter; id++ ) {
    TGeoNode* daughternode = volume->GetNode(id);
    cout << preface.c_str() << "  " << id << ")" << daughternode->GetName() 
         << endl;
    cout << preface.c_str() << "  ";
    daughternode->GetMatrix()->Print();
    TGeoVolume* daughtervol = daughternode->GetVolume();
    this -> DumpVolume(daughtervol->GetName(),preface+"  ");
  }
  
}

//_____________________________________________________________________________
std::string GeoGeometry::GetGeoCompatibleName(std::string name) {
  // Public Static method to return acceptable node or volume name
  // Symbol "/" is replaced with "f" to avoid conflict with "/" used in paths
  // Symbol "-" is replaced with "h" to avoid conflict with "-" used in
  // composite shapes. 

  std::string okayname = name;
  // Name cannot include "/" because of ambiguity with "/" used in
  // path to nodes. Replace "/" with "f".
  unsigned int ipos = okayname.find("/");
  while ( ipos != std::string::npos ) {
    okayname.replace(ipos,1,"f");
    ipos = okayname.find("/");
  }

  // Name cannot include "-" because of ambiguity with "-" used in
  // composite shapes. Replace "-" with "h".
  ipos = okayname.find("-");
  while ( ipos != std::string::npos ) {
    okayname.replace(ipos,1,"h");
    ipos = okayname.find("-");
  }
  
  return okayname; 
}

//_____________________________________________________________________________
TVector3 GeoGeometry::GetHallExtentMin() const {
  // Purpose:: Return the min (x,y,z) of the detector hall in global
  // (MARS) coordinates.

  UpdateGlobalManager();
  
  gGeoManager -> cd(fHallPath.c_str());
  TGeoNode* hallNode = gGeoManager->GetCurrentNode();
  TGeoVolume* hallVol = hallNode->GetVolume();

  TGeoBBox* hallBox = dynamic_cast<TGeoBBox*>(hallVol->GetShape());
  
  Double_t lxyz[3] = {-(hallBox->GetDX()),-(hallBox->GetDY()),
                      -(hallBox->GetDZ())};
  Double_t gxyz[3] = {0};

  gGeoManager->LocalToMaster(lxyz,gxyz);
  
  return TVector3(gxyz[0],gxyz[1],gxyz[2]);

}

//_____________________________________________________________________________
TVector3 GeoGeometry::GetHallExtentMax() const {
  // Purpose:: Return the max (x,y,z) of the detector hall in global
  // (MARS) coordinates.

  UpdateGlobalManager();
  
  gGeoManager -> cd(fHallPath.c_str());
  TGeoNode* hallNode = gGeoManager->GetCurrentNode();
  TGeoVolume* hallVol = hallNode->GetVolume();

  TGeoBBox* hallBox = dynamic_cast<TGeoBBox*>(hallVol->GetShape());
  
  Double_t lxyz[3] = {+(hallBox->GetDX()),+(hallBox->GetDY()),
                      +(hallBox->GetDZ())};
  Double_t gxyz[3] = {0};
  
  gGeoManager->LocalToMaster(lxyz,gxyz);
  
  return TVector3(gxyz[0],gxyz[1],gxyz[2]);

}

//_____________________________________________________________________________
GeoSteelPlnNode* GeoGeometry::GetNearestSteelPlnNode(Double_t z) const {
  // Public. Purpose:: Retrieve the nearest steel node for a given z

  UpdateGlobalManager();

  MSG("Geo",Msg::kDebug) << "GeoGeometry::GetNearestSteelPlnNode to z " 
                         << z << "." << endl;
  assert(fGeoManager);

  const std::map<Float_t,GeoSteelPlnNode*>& stnodemap = GetSteelPlnNodeMap();

  // Steel planes in map are all in main detector block
  // Ordering of planes is lo to hi
  std::map<Float_t,GeoSteelPlnNode*>::const_iterator steelplnItr;
  
  // Lower bound returns steel plane with z value greater than or equal to key
  steelplnItr = stnodemap.lower_bound(z);
  GeoSteelPlnNode* loplnNode = 0;
  Double_t         loplnz0   = -9999.;
  GeoSteelPlnNode* hiplnNode = 0;
  Double_t         hiplnz0   = -9999.;
  if ( steelplnItr != stnodemap.end() ) {
    hiplnNode = steelplnItr->second;
    hiplnz0   = steelplnItr->first;
  }
  if ( steelplnItr != stnodemap.begin() ) {
    steelplnItr--;
    loplnNode = steelplnItr->second;
    loplnz0   = steelplnItr->first;
  }
  if ( hiplnNode == 0 ) return loplnNode;
  else if ( loplnNode == 0 ) return hiplnNode;
  
  // Bracketed by two planes, determine which is closer
  // avoid additional (slow) coordinate transformations by using
  // information stored in the map
  return ( TMath::Abs(z-loplnz0) < TMath::Abs(z-hiplnz0) ) ? loplnNode : hiplnNode;


}

//_____________________________________________________________________________
PlexPlaneId GeoGeometry::GetPlaneIdFromZ(Double_t z) const {
  // Purpose:: Retrieve detector plane corresponding to the first
  // plane in the main detector block with a back face (high-z side)
  // "downstream" of z (greater than z).  The exception is that if z 
  // is greater than back face of highest-z plane in detector, return 
  // plexplaneid of last plane, even though this plane's high-z face is 
  // not downstream of z.  This behavior is that used in UgliGeometry. 
  // Use PlexPlaneId::IsValid() to determine validity of returned PlexPlaneId.
  // Will only return invalid plexplaneid if no planes in main detector
  // block.  
 
  UpdateGlobalManager();
  
  MSG("Geo",Msg::kDebug) << "GeoGeometry::GetPlaneIdFromZ " << z << "."<< endl;
  assert(fGeoManager);

  PlexPlaneId defplex;  // default plex is invalid
  
  PlaneMapConstItr plnItr;
  // Move through the plane map from low to high z
  Double_t xyz_global[3] = {0,0,z};
  Double_t xyz_local[3] = {0};
      
  // fPlaneMap has scint & steel planes 
  for ( plnItr = fPlaneMap.begin();  plnItr != fPlaneMap.end(); ++plnItr) {
    if ( (plnItr->first).IsVetoShield() ) continue; // No shield
    GeoPlnNode* plnNode = dynamic_cast<GeoPlnNode*>(plnItr -> second);
    plnNode->GlobalToLocal(xyz_global,xyz_local);
    Float_t dz  = dynamic_cast<TGeoBBox*>(plnNode->GetVolume()->GetShape())
                  ->GetDZ();
    if ( xyz_local[2] < dz ) {
      return plnItr->first;
    }
    defplex = plnItr->first; // default plex is highest-z pln in main block
  }

  return defplex;

}

//_____________________________________________________________________________
const vector<GeoPlnNode*>& GeoGeometry::GetPlnNodePtrVector() const {
  // Public. Purpose:: Return collection of ptrs for all pln nodes in detector

  MSG("Geo",Msg::kDebug) << "GetPlnNodePtrVector" << endl;

  UpdateGlobalManager();

  if ( fPlnNodes.empty() ) {
    PlaneMapConstItr plnItr;
    for ( plnItr = fPlaneMap.begin(); plnItr != fPlaneMap.end(); ++plnItr) {
      GeoPlnNode* plnNode = plnItr->second;
      fPlnNodes.push_back(plnNode);
    }
  }
  
  return fPlnNodes;

}

//_____________________________________________________________________________
GeoScintPlnNode* GeoGeometry::GetScintPlnNode(PlexPlaneId planeid) const {
  // Public. Purpose:: Retrieve the node for a particular plane

  MSG("Geo",Msg::kDebug) << "GetScintPlnNode for plane " << planeid.AsString() 
                         << endl;

  UpdateGlobalManager();

  GeoScintPlnNode* plnnode = 0;
  if ( planeid.IsSteel() ) {
    MSG("Geo",Msg::kWarning) 
      << "GetScintPlnNode called with steel plane id " 
      << planeid.AsString() << endl;
    return 0;
  }
  
  PlaneMapConstItr citr = fPlaneMap.find(planeid);
  if ( citr != fPlaneMap.end() ) {
    plnnode = dynamic_cast<GeoScintPlnNode*>(citr->second);
  }
  
  return plnnode;

}

//_____________________________________________________________________________
const vector<GeoScintPlnNode*>& GeoGeometry::GetScintPlnNodePtrVector() const {
  // Purpose:: Return collection of ptrs for all scint pln nodes in detector

  MSG("Geo",Msg::kDebug) << "GetScintPlnNodePtrVector" << endl;

  UpdateGlobalManager();

  if ( fScintPlnNodes.empty() ) {
    PlaneMapConstItr plnItr;
    for ( plnItr = fPlaneMap.begin(); plnItr != fPlaneMap.end(); ++plnItr ){
      if ( !(plnItr->first).IsSteel() ) {
        fScintPlnNodes.push_back(
                              dynamic_cast<GeoScintPlnNode*>(plnItr->second));
      }
    }
  }

  return fScintPlnNodes;

}

//_____________________________________________________________________________
GeoSteelPlnNode* GeoGeometry::GetSteelPlnNode(PlexPlaneId planeid) const {
  // Public. Purpose:: Retrieve the node for a particular plane

  MSG("Geo",Msg::kDebug) << "GetSteelPlnNode for plane " << planeid.AsString() 
                         << endl;

  UpdateGlobalManager();

  GeoSteelPlnNode* plnnode = 0;
  if ( !planeid.IsSteel() ) {
    MSG("Geo",Msg::kWarning) 
      << "GetSteelPlnNode called with non steel plane id " 
      << planeid.AsString() << endl;
    return 0;
  }
  
  PlaneMapConstItr citr = fPlaneMap.find(planeid);
  if ( citr != fPlaneMap.end() ) {
    plnnode = dynamic_cast<GeoSteelPlnNode*>(citr->second);
  }
  
  return plnnode;

}

//_____________________________________________________________________________
const vector<GeoSteelPlnNode*>& GeoGeometry::GetSteelPlnNodePtrVector() const {
  // Public. Purpose:: Return collection of ptrs for all steel pln nodes 
  // in detector

  MSG("Geo",Msg::kDebug) << "GetSteelPlnNodePtrVector" << endl;

  UpdateGlobalManager();

  if ( fSteelPlnNodes.empty() ) {
    PlaneMapConstItr plnItr;
    for ( plnItr = fPlaneMap.begin(); plnItr != fPlaneMap.end(); ++plnItr ){
      if ( (plnItr->first).IsSteel() ) {
        fSteelPlnNodes.push_back(
                              dynamic_cast<GeoSteelPlnNode*>(plnItr->second));
      }
    }
  }
  
  return fSteelPlnNodes;

}

//_____________________________________________________________________________
const std::map<Float_t,GeoSteelPlnNode*>& GeoGeometry::GetSteelPlnNodeMap() 
                                                                     const {
  // Public. Purpose:: Return collection of ptrs for all steel pln nodes 
  // in main body of detector, key'ed by z-position. 

  MSG("Geo",Msg::kDebug) << "GetSteelPlnNodeMap" << endl;

  UpdateGlobalManager();
  
  if ( fSteelPlnNodeMap.empty() ) {
    PlaneMapConstItr plnItr;
    for ( plnItr = fPlaneMap.begin(); plnItr != fPlaneMap.end(); plnItr++ ){
      if ( (plnItr->first).IsSteel() && !(plnItr->first).IsVetoShield()) {
        fSteelPlnNodeMap.insert(make_pair(plnItr->second->GetZ0(),
                             dynamic_cast<GeoSteelPlnNode*>(plnItr->second)));
      }
    }
  }
  
  return fSteelPlnNodeMap;

}

//_____________________________________________________________________________
GeoStripNode* GeoGeometry::GetStripNode(const PlexStripEndId& seid) const {
  // Public. Purpose:: Retrieve the node for a particular strip
 
  MSG("Geo",Msg::kDebug) << "GetStripNode for strip "  << seid.AsString() 
                         << endl;

  UpdateGlobalManager();

  GeoStripNode* stripnode = 0;

  // Retrieve the plane node and then the strip node
  GeoScintPlnNode* planenode = this -> GetScintPlnNode(seid);
  if ( planenode ) {
    stripnode = planenode -> GetStripNode(seid);
  }

  return stripnode;

}

//_____________________________________________________________________________
void GeoGeometry::GetTransverseExtent(PlaneView::PlaneView_t view,
                                      Double_t& tmin, Double_t& tmax) const {
  // Public. Purpose: Retrieve transverse extent of detector
  // based on detector type and form.  These values are currently hardwired
  // as in UgliGeometry::GetTransverseExtent

  MSG("Geo",Msg::kDebug) << "GetTransverseExtent for view "  
                         << PlaneView::AsString(view) << endl;

  UpdateGlobalManager();

  tmin = -999; tmax=-999;

  Detector::Detector_t detector 
    = (Detector::Detector_t)fVldRange.GetDetectorMask();
  
  switch (detector) {
  case Detector::kNear:
    switch (view ) {
    case PlaneView::kU:
      tmin = -2.10 * fScale;
      tmax = +2.75 * fScale;
      break;
    case PlaneView::kV:
      tmin = -2.75 * fScale;
      tmax = +2.10 * fScale;
      break;
    default:
      MSG("Geo",Msg::kError)
        << "GeoGeometry::GetTransverseExtent undefined for "
        << Detector::AsString(detector) << " view "
        << PlaneView::AsString(view) << endl;
      break;
    }
  case Detector::kFar:
    tmin = -4.0 * fScale;
    tmax = +4.0 * fScale;
    break;
  case Detector::kCalib:
    tmin = -0.5 * fScale;
    tmax = +0.5 * fScale;
    break;
  default:
    MSG("Geo",Msg::kError) 
      << "GeoGeometry::GetTransverseExtent undefined for "
      << Detector::AsString(detector) << endl;
    break;
  }

  Float_t fuzz = 0.5*Geo::kStripWidth*fScale;
  tmin = tmin - fuzz;
  tmax = tmax + fuzz;
  
  return;
    
}

//_____________________________________________________________________________
void GeoGeometry::GetTransverseExtent(PlaneView::PlaneView_t view,
                                      Float_t& tmin, Float_t& tmax) const {
  // Public. Purpose: Retrieve transverse extent of detector in global (MARS)
  // coordinates.

  UpdateGlobalManager();
  
  Double_t tmind,tmaxd;
  GetTransverseExtent(view,tmind,tmaxd);
  tmin = (Float_t)tmind;
  tmax = (Float_t)tmaxd;

  return;
    
}

//_____________________________________________________________________________
void GeoGeometry::GetZExtent(Double_t& zmin, Double_t& zmax, Int_t isup) const{
  // Public. Purpose: Retrieve zextent of specified isup supermodule in MARS
  // coordinates. If isup = -1 (default), returns zextent of entire detector.
 
  MSG("Geo",Msg::kDebug) << "GetZExtent for isup "  << isup << endl;

  UpdateGlobalManager();
 
  // The boundaries are calculated during geometry build 
  zmax =-999999;
  zmin =+999999; 
  if ( isup < 0 ) {
    zmin = fSMZMin[0];
    if ( fSMZMax[1] > 0 ) zmax = fSMZMax[1];
    else zmax = fSMZMax[0];
  }
  else if ( isup < 2 ) {
    zmin = fSMZMin[isup];
    zmax = fSMZMax[isup];
  }
  else {
    MSG("Geo",Msg::kWarning) << "GetZExtent failed. SuperModule "
                             << isup << " undefined." << endl;
  }

  return;
  
   
}

//_____________________________________________________________________________
void GeoGeometry::GetZExtent(Float_t& zmin, Float_t& zmax, Int_t isup) const {
  // Public. Purpose:: Retrieve zextent of specified isup supermodule. If 
  // isup = -1 (default), returns zextent of entire detector.

  UpdateGlobalManager();
  
  Double_t zmind,zmaxd;
  GetZExtent(zmind,zmaxd,isup);
  zmin = (Float_t)zmind;
  zmax = (Float_t)zmaxd;

  return;
}

//_____________________________________________________________________________
TGeoMedium* GeoGeometry::GetMedium(const char* name) const {
  // Public. Purpose:: Return ptr to medium of given name.  If geometry
  // is not closed, such that this method is being called during the
  // geometry build stage, will abort with fatal message if no medium
  // of given name is found.  Otherwise, will return null ptr if
  // medium not found.
  //  

  UpdateGlobalManager();

  TGeoMedium* medium = fGeoManager->GetMedium(name);
  if ( !medium && !(fGeoManager->IsClosed()) ) {
    MSG("Geo",Msg::kFatal) << "GeoGeometry::GetMedium did not find medium\n"
                           << name << " during geometry build stage. Abort."
                           << endl;
    abort();
  }
  
  return medium;

}

//_____________________________________________________________________________
bool GeoGeometry::IsCompatible(const VldContext& vldc, 
                               Geo::EAppType apptype) const {
  // Public. Test if vldc and apptype are in range of this geometry

  UpdateGlobalManager();

  bool iscompatible = false;
  if ( fVldRange.IsCompatible(vldc) && fAppType == apptype ) 
                                             iscompatible = true;
  return iscompatible; 
}

//_____________________________________________________________________________
void GeoGeometry::ls(Option_t* option) const {
  // Public. list components of this geometry
  // Argument: option string is a concatenated list of char options, where:
  //           "m" => print list materials
  //           "M" => print list of media.  Default won't print swim media
  //                  Override with configopt == "a", e.g. ls("Ma").
  //           "d" => print MediumMap detector components and associated 
  //                  mediums and swim methods 
  //           "v" => print list of volumes 
  //           "a" => configuration option "all".  If present, may toggle
  //                  on a greater dump of information.
  // Example: An option argument to specify listing of materials and mediums
  //          is "mM"

  UpdateGlobalManager();
  PrintHeader();

  TString opt = option;
  bool isAll = false;
  if ( opt.Contains("a") ) isAll = true;
  
  if (opt.Contains("m")) {
    // Can't use ls() for TGeoMaterial because not supported
    fGeoManager->GetListOfMaterials()->Print(); 
  }
  if (opt.Contains("M")) {
    // Can't use ls() or Print() for TGeoMedium because not supported
    GeoMedium::PrintHeader();  // print media table header
    TIter medItr(fGeoManager->GetListOfMedia());
    GeoMedium* medium = 0;  
    while ( ( medium = dynamic_cast<GeoMedium*>(medItr()) ) ) {
      if ( !isAll ) {
        std::string medName = medium->GetName();
        if ( medName.find("_SWIM") != std::string::npos ) continue;
      }
      medium -> Print();
    }
  }
  if (opt.Contains("d")) {
    // Print list of detector components and associated mediums and
    // swim methods
    GetMediumMap().Print();
  }
  if (opt.Contains("v")) {
    // Print list of volumes in this geometry
    TIter volItr(fGeoManager->GetListOfUVolumes());
    TGeoVolume* volume = 0;
    while ( ( volume = dynamic_cast<TGeoVolume*>(volItr()) ) ) {
      cout << setiosflags(ios::left) << setw(30) << volume->GetName() 
           << setw(30) << volume->GetMedium()->GetName()
           << endl;
    }
  }    
  
}

//_____________________________________________________________________________
void GeoGeometry::Print(Option_t* option) const {
   // Public. Print something about this geometry (name + vldrange)

  UpdateGlobalManager();

  std::string stoption(option);
  if (stoption != "") DumpVolume(stoption);
  PrintHeader();
  
}

//_____________________________________________________________________________
void GeoGeometry::PrintHeader(Option_t* /* option */) const {
  // Print header info (name + vldrange) about this geometry.

  cout << "GeoGeometry " << fGeoManager->GetName() << " " 
       << fGeoManager->GetTitle() << endl;
  cout << "Has " << CountRef() << " reference" << ((CountRef()==1)?"":"s")
       << " VldRange:" << fVldRange.AsString("acs-") << endl;
  
}

//_____________________________________________________________________________
void GeoGeometry::SwitchMedia(bool toswim) {
  // Method to geometry to/from swim media to default media
  // Swim media are used during GeoSwimmer particle transport and
  // have a different set of physics process flags/energy threshold
  // cuts assigned to them when used with a concrete VMC.

  Int_t ntoadd = 0;
  Int_t nmedia = fGeoManager->GetListOfMedia()->GetSize();
  if ( toswim ) {
    if ( fIsSwimMedia ) return;  // already set up for swim media
    ntoadd = nmedia/2 - 1; // number to add to current medId to get swim twin  
  }
  else {
    if ( !fIsSwimMedia ) return; // already set up for default media
    ntoadd = -nmedia/2 - 1; // number to add to current medId to get def twin
  }
  
  TIter nextvol(fGeoManager->GetListOfVolumes());
  TGeoVolume* vol = 0;
  while ( ( vol = (TGeoVolume*)nextvol() ) ) {
    TGeoMedium* currentMed = vol -> GetMedium();
    Int_t twinMedId = currentMed->GetId() + ntoadd;
    TGeoMedium* twinMed = (TGeoMedium*)
                     (fGeoManager->GetListOfMedia()->At(twinMedId));
    vol -> SetMedium(twinMed);
  }

  fIsSwimMedia = toswim;

}

//_____________________________________________________________________________
void GeoGeometry::BuildAll(const VldContext &vldc) {
  // Private.  Build materials & geometry

  UpdateGlobalManager();

  MSG("Geo",Msg::kDebug) << "GeoGeometry BuildAll " << endl;
  assert(fGeoManager);

  // Get the VldRange that fits this context
  BuildVldRange(vldc);

  // Build commonly used rotation matrices
  BuildRotations();
  
  // Build Detector Geometry
  BuildGeometry(vldc);

}

//_____________________________________________________________________________
void GeoGeometry::BuildVldRange(const VldContext &vldc) {
  // Private. Build VldRange from VldContext

  UpdateGlobalManager();

  MSG("Geo",Msg::kDebug) << "GeoGeometry::BuildVldRange for VldContext " 
                         << vldc << endl;
  
  Detector::Detector_t det = vldc.GetDetector();
  bool isunknown = (Detector::kUnknown == det);
  const char* src = "DBI";
  Int_t detmask = Detector::FullMask();
  Int_t simmask = SimFlag::FullMask();
  if ( isunknown ) { 
    src = "Fake (unknown detector)";
    detmask = det;
    simmask = vldc.GetSimFlag();
  }
  
  // Start the VldRange covering all of time and then trim it down
  VldTimeStamp starttime = VldTimeStamp::GetBOT();
  VldTimeStamp endtime = VldTimeStamp::GetEOT();
  fVldRange = VldRange(detmask,simmask,starttime,endtime,src);

  // special return for unknown detector
  if (isunknown) return;

  DbiResultPtr<UgliDbiSteelPln> steelTbl(vldc);
  TrimVldRange("UgliDbiSteelPln",steelTbl.GetValidityRec());

  DbiResultPtr<UgliDbiScintPlnStruct> scintStructTbl(vldc);
  TrimVldRange("UgliDbiScintPlnStruct",scintStructTbl.GetValidityRec());

  DbiResultPtr<UgliDbiScintPln> scintTbl(vldc);
  TrimVldRange("UgliDbiScintPln",scintTbl.GetValidityRec());

  DbiResultPtr<UgliDbiScintMdlStruct> mdlStructTbl(vldc);
  TrimVldRange("UgliDbiScintMdlStruct",mdlStructTbl.GetValidityRec());

  DbiResultPtr<UgliDbiScintMdl> mdlTbl(vldc);
  TrimVldRange("UgliDbiScintMdl",mdlTbl.GetValidityRec());

  DbiResultPtr<UgliDbiStripStruct> stripStructTbl(vldc);
  TrimVldRange("UgliDbiStripStruct",stripStructTbl.GetValidityRec());

  DbiResultPtr<UgliDbiStrip> stripTbl(vldc);
  TrimVldRange("UgliDbiStrip",stripTbl.GetValidityRec());
  
  MSG("Geo",Msg::kDebug) << "GeoGeometry::BuildVldRange built VldRange " 
                         << fVldRange << endl;
  
}

//_____________________________________________________________________________
void GeoGeometry::TrimVldRange(const char* tblName,
                               const DbiValidityRec* dbivrec) {
  // Private. Trim VldRange based on DbiValidityRec range info

  UpdateGlobalManager();

  if ( dbivrec ) {
    fVldRange.TrimTo(dbivrec->GetVldRange());
  }
  else {
    MSG("Geo",Msg::kWarning) << "No DbiValidityRec for table " << tblName 
                             << endl;
  }
  
}
   
//_____________________________________________________________________________
void GeoGeometry::BuildRotations() {
  // Build commonly used rotation matrices

  // Rotation about the z-axis by +45 degrees
  fRot45 = new TGeoRotation("rot45",90,45,90,135,0,0);
  fRot45 -> RegisterYourself(); // do this so that it can be found later


  // Rotation about the z-axis by -45 degrees
  fRot315 = new TGeoRotation("rot315",90,315,90,45,0,0);
  fRot315 -> RegisterYourself(); // do this so that it can be found later

}


//_____________________________________________________________________________
void GeoGeometry::BuildGeometry(const VldContext& vldc) {
  // Private. Build this geometry's hierarchy of nodes

  UpdateGlobalManager();
  
  MSG("Geo",Msg::kDebug) << "GeoGeometry::BuildGeometry" << endl;
  assert(fGeoManager);

  const GeoMediumMap& medMap = GetMediumMap();
  
  if (vldc.GetDetector() == Detector::kUnknown) {
    // build a minimalistic geometry
    TGeoVolume* volMARS = fGeoManager -> MakeBox("MARS",
                                   medMap.GetMedium(Geo::kMars),1.2,1.2,1.2);
    TGeoVolume* volLINR = fGeoManager -> MakeBox("LINR",
                                   medMap.GetMedium(Geo::kLinr),1.1,1.1,1.1);
    volMARS -> AddNode(volLINR,1,gGeoIdentity);
    TGeoVolume* volHALL = fGeoManager -> MakeBox("HALL",
                                   medMap.GetMedium(Geo::kHall),1.0,1.0,1.0);
    volLINR -> AddNode(volHALL,1,gGeoIdentity);
    fGeoManager -> SetTopVolume(volMARS);
    return;
  }

  const UgliDbiTables& ugliDbiTables = UgliDbiTables(vldc);
  const UgliDbiGeometry* ugliDbiGeo = ugliDbiTables.GetDbiGeometry();
  if (!ugliDbiGeo) {
    MSG("Geo",Msg::kFatal) 
       << "No UgliDbiGeometry for " << vldc << endl;
    abort();
  }

  // Start with the hall, this is the widest view stored in db 
  const Float_t* hmin = ugliDbiGeo->GetHallXYZMin();
  const Float_t* hmax = ugliDbiGeo->GetHallXYZMax();

  // For far detector, pad x-dimension because otherwise shield is cut-off
  Double_t xpad = 0.*fScale;
  if ( vldc.GetDetector() == Detector::kFar ) xpad = 0.02*fScale;
  // For near detector, pad y-dim by 1 cm because otherwise scint plane
  // polygon extrudes floor in imperfect geometry. Fix me.
  Double_t ypad = 0.*fScale;
  if ( vldc.GetDetector() == Detector::kNear ) ypad = 0.01*fScale;
  
  Double_t hminsc[3] = {hmin[0]*fScale,hmin[1]*fScale,hmin[2]*fScale};
  Double_t hmaxsc[3] = {hmax[0]*fScale,hmax[1]*fScale,hmax[2]*fScale};

  // x/y/zhsiz is half-width in each dimension of hall
  Double_t xhsiz = 0.5*(hmaxsc[0] - hminsc[0]) + xpad;
  Double_t yhsiz = 0.5*(hmaxsc[1] - hminsc[1]) + ypad;
  Double_t zhsiz = 0.5*(hmaxsc[2] - hminsc[2]);

  // x/y/z0hall is center of hall
  Double_t x0hall = 0.5*(hmaxsc[0] + hminsc[0]);
  Double_t y0hall = 0.5*(hmaxsc[1] + hminsc[1]);
  Double_t z0hall = 0.5*(hmaxsc[2] + hminsc[2]);

  MSG("Geo",Msg::kDebug) << "Hall halfwidths " << xhsiz << ","
                         << yhsiz << "," << zhsiz << " placed at "
                         << x0hall << "," << y0hall << "," << z0hall << endl;
  

  // Hall
  TGeoVolume* volHALL   = fGeoManager->MakeBox("HALL",
                             medMap.GetMedium(Geo::kHall),xhsiz,yhsiz,zhsiz);
  volHALL-> SetVisibility(kTRUE);
  volHALL-> SetLineColor(kGreen);

  // Liner pad is box with halfwidth dimensions equal
  // to the hall halfwidths + the liner thickness (1 m)
  Detector::Detector_t dettype = (Detector::Detector_t)vldc.GetDetector();
  const Double_t linerpad = 1.0*fScale; // thickness concrete liner (1 m)
  
  TGeoVolume* volLINR     = fGeoManager->MakeBox("LINR",
                                         medMap.GetMedium(Geo::kLinr),
                                         xhsiz+linerpad,yhsiz+linerpad,
                                         zhsiz+linerpad);
  volLINR-> SetVisibility(kTRUE);
  volLINR-> SetLineColor(kYellow);

  // Build MARS appropriate to detector type.
  TGeoVolume* volMARS = 0;
  switch ( dettype ) {
  case Detector::kFar:
    volMARS = BuildFarMARS();
    break;
  case Detector::kNear:
    volMARS = BuildNearMARS();
    break;
  default:
    // Rock shell default is
    // rock padding, minimum 5 m thick onto half hall dimensions
    const Double_t rockpad = 5.0*fScale; // thickness of considered rock (5 m)
    Double_t absxmax = TMath::Max(TMath::Abs(hminsc[0]-xpad),
                                  TMath::Abs(hmaxsc[0]+xpad));
    Double_t absymax = TMath::Max(TMath::Abs(hminsc[1]-ypad),
                                TMath::Abs(hmaxsc[1]+ypad));
    Double_t abszmax = TMath::Max(TMath::Abs(hminsc[2]),
                                TMath::Abs(hmaxsc[2]));
    volMARS = fGeoManager->MakeBox("MARS",medMap.GetMedium(Geo::kMars),
                                    absxmax+rockpad,
                                    absymax+rockpad,abszmax+rockpad);
    break;
  }

  volMARS -> SetVisibility(kTRUE);
  volMARS -> SetLineColor(kBlack); 

  // The top level volume MARS will be assigned node path "/MARS_1" by root
  fGeoManager -> SetTopVolume(volMARS);
  // GeoNodes are owned by the TGeoManager
  std::string nodename = "LINR_1";
  std::string globalpath = "/MARS_1/" + nodename;
  volMARS -> AddNode(volLINR,1,
                     new TGeoTranslation("LINR",x0hall,y0hall,z0hall));
   
  nodename = "HALL_1";
  globalpath += "/" + nodename;
  fHallPath = globalpath;
  volLINR -> AddNode(volHALL,1,gGeoIdentity);

  BuildDetector(vldc,volHALL);
  
}

//_____________________________________________________________________________
void GeoGeometry::BuildDetector(const VldContext& vldc,
                                TGeoVolume* hallVol) {
   // Private. Build this geometry's steel+scint planes + coil

  UpdateGlobalManager();

  MSG("Geo",Msg::kDebug) << "GeoGeometry::BuildDetector " << endl;
  assert(fGeoManager);

  DbiResultPtr<UgliDbiSteelPln> steelTbl(vldc);

  BuildPlanePairVolumes(vldc); // pre-build plane pair volumes

  // Build shield group volumes and install them as nodes in hall
  // and install shield plane pair nodes in group volumes
  if ( vldc.GetDetector() == Detector::kFar ) 
                             fGeoShield.BuildGroupNodes(hallVol);
  
  // Plane db positions are stored relative to mars so need
  // hall global position to recalculate
  // Translation taken from liner, since hall placed with identity
  // matrix relative to liner
  TGeoNode* linrNode = fGeoManager->GetVolume("MARS")->GetNode("LINR_1");
  const Double_t* h0 = linrNode -> GetMatrix() -> GetTranslation();
  
  Int_t plnMin[2] = {-1};
  Int_t plnMax[2] = {-1};
  Int_t nSM = GetSMPlaneLimits(vldc,plnMin,plnMax);

  std::string minPlaneName[2];
  std::string maxPlaneName[2];
  Int_t minPlane[2] = {9999,9999};
  Int_t maxPlane[2] = {-9999,-9999};
  
  // order construction by plane #, STL set will sort them
  std::set<PlexPlaneId> steelset;
  for (unsigned int irow=0; irow < steelTbl.GetNumRows(); ++irow) {
    const UgliDbiSteelPln* stRow = steelTbl.GetRow(irow);
    PlexPlaneId steelId = stRow -> GetPlaneId();
    steelset.insert(steelId);
  }
  
  std::set<PlexPlaneId>::const_iterator steelitr = steelset.begin();
  for ( ; steelitr != steelset.end() ; steelitr++) {
    PlexPlaneId steelId = *steelitr;
    const UgliDbiSteelPln* stRow = steelTbl.GetRowByIndex(steelId.GetPlane());
    if ( (vldc.GetDetector() ==  Detector::kFar && !steelId.IsVetoShield()) ) {
      // Dispatch some anomalies in the database
       if ( steelId.GetPlane() > 485 ) continue; // because db has 0-497
       // Skip planes that have intentionally been placed at unphysical
       // locations
       if ( TMath::Abs(stRow->GetX0() ) > 50. ) continue;
    }

    // This is used just for reporting info at end of build
    Int_t vsId = 0;
    std::string stName = GetGeoCompatibleName(steelId.AsString());
    if ( steelId.IsVetoShield() ) vsId = 1;
    if ( steelId.GetPlane() < minPlane[vsId] ) {
        minPlane[vsId] = steelId.GetPlane();
        minPlaneName[vsId] = stName;
    }
    if ( steelId.GetPlane() > maxPlane[vsId] ) {
        maxPlane[vsId] = steelId.GetPlane();
        maxPlaneName[vsId] = stName;
    }
    
    // Veto shield built separately
    if ( vldc.GetDetector() == Detector::kFar && steelId.IsVetoShield() ) 
                                                              continue;
    
    Double_t offxyz[3] = {h0[0],h0[1],h0[2]};
       
    // The pair volumes have been pre-built
    std::string pairName = GetGeoCompatibleName(steelId.AsString("b"));
    TGeoVolume* pairVol = fGeoManager -> GetVolume(pairName.c_str());
    if ( !pairVol ) {
      MSG("Geo",Msg::kFatal) 
         << "Unable to find pre-built pairvol for volume "
         << pairName.c_str() << endl;
      abort();
    }
    
    // Position pair bounding box 
    TGeoRotation* stRotMatrix = new TGeoRotation(pairName.c_str(),
                                stRow->GetThetaDeg(0),stRow->GetPhiDeg(0),
                                stRow->GetThetaDeg(1),stRow->GetPhiDeg(1),
                                stRow->GetThetaDeg(2),stRow->GetPhiDeg(2));
    stRotMatrix -> RegisterYourself();
    
    // Plane DB translation is relative to MARS, move to SM or hall local 
    Float_t x0pair = (stRow->GetX0())*fScale - offxyz[0];
    Float_t y0pair = (stRow->GetY0())*fScale - offxyz[1];

    // If plane is first plane in supermodule than thickness of pair
    // box is just "Thickness".  If it's part of vetoshield, than 
    // thickness of pair box is just "TotalZ", else pair box thickness 
    // is set by distance from back face of
    // of current steel plane to back face of previous steel plane. z0pair
    // is defined as the position of the center of this pair box in z.
    Int_t planeno = steelId.GetPlane();
    Float_t z0pair = 0;
    if ( vldc.GetDetector() == Detector::kCalDet || steelId.IsVetoShield() ) 
      z0pair =(stRow->GetZBack()-0.5*(stRow->GetTotalZ()))*fScale - offxyz[2];
    else
      z0pair =(stRow->GetZBack()-0.5*(stRow->GetThickness()))*fScale-offxyz[2];
          
    if ( !steelId.IsVetoShield() ) {
      if ( planeno != plnMin[0] ) {
        if ( nSM == 1 || (nSM > 1 && planeno != plnMin[1]) ) {
          const UgliDbiSteelPln* stPrevRow = steelTbl.GetRowByIndex(planeno-1);
          z0pair = (stRow->GetZBack()-0.5*(stRow->GetZBack()
                                  -stPrevRow->GetZBack()))*fScale - offxyz[2];
        }
      }
    }
    
    TGeoCombiTrans* pairMatrix = new TGeoCombiTrans(pairName.c_str(),
                                    x0pair,y0pair,z0pair,stRotMatrix);
    MSG("Geo",Msg::kDebug) << "Placing plane pair " << pairName.c_str()
                           << " at " << x0pair << ","
                           << y0pair << "," << z0pair 
                           << "\n rotmatrix:(thetax/y/z,phix/y/z) (" 
                           << stRow->GetThetaDeg(0) << ","  
                           << stRow->GetPhiDeg(0) << "),("
                           << stRow->GetThetaDeg(1) << ","  
                           << stRow->GetPhiDeg(1) << "),("
                           << stRow->GetThetaDeg(2) << ","  
                           << stRow->GetPhiDeg(2) << ")"
                           << endl;

    // Place plane pairs directly in hall
    hallVol -> AddNode(pairVol,1,pairMatrix);
  } // end of loop over all plane db rows


  MSG("Geo",Msg::kInfo) << "GeoGeometry built detector planes "
                        << minPlaneName[0].c_str() << " to " 
                        << maxPlaneName[0].c_str() << "." << endl;
  if ( !minPlaneName[1].empty() ) {
    MSG("Geo",Msg::kInfo) << "GeoGeometry built shield planes "
                          << minPlaneName[1].c_str() << " to " 
                          << maxPlaneName[1].c_str() << "." << endl;
  }
  
  // Build outer coil after all planes are positioned so that detector
  // limits are known
  Detector::Detector_t dettype = (Detector::Detector_t)vldc.GetDetector();
  switch ( dettype ) {
  case Detector::kFar:
    BuildFarCoil(hallVol);
    break;
  case Detector::kNear:
    BuildNearCoil(hallVol);
    break;
  default:
    break;
  }

  return;
  
}

//_____________________________________________________________________________
void GeoGeometry::BuildPlanePairVolumes(const VldContext& vldc) {
   // Private. Build this geometry's steel+scint plane pair volumes

  UpdateGlobalManager();

  MSG("Geo",Msg::kDebug) << "GeoGeometry::BuildPlanePairVolumes " << endl;
  assert(fGeoManager);

  const GeoMediumMap& medMap = GetMediumMap();

  DbiResultPtr<UgliDbiSteelPln> steelTbl(vldc);
  const UgliDbiTables& ugliDbiTables = UgliDbiTables(vldc);

  // Pre-build strip volumes
  BuildStripVolumes(vldc);

  Int_t plnMin[2] = {-1};
  Int_t plnMax[2] = {-1};
  Int_t nSM = GetSMPlaneLimits(vldc,plnMin,plnMax);
  
  // Build coil segment volume to be used for air gaps between scint & steel
  // planes.  This is built as TGeoVolumeMulti, and can be built once and
  // used to fill every air gap in detector.
  TGeoVolume* coilVol = BuildCoilAirGapVolume(vldc); 
  TGeoRotation* rot45 = GetRot45();
  Detector::Detector_t dettype = (Detector::Detector_t)vldc.GetDetector();
  
  for (unsigned int irow=0; irow < steelTbl.GetNumRows(); ++irow) {

    // Build steel and scint plane volumes first, then position so
    // that we know pair volume size.

    // Steel plane
    GeoSteelPlnVolume* stplnVolume = 0;
    const UgliDbiSteelPln* stRow = steelTbl.GetRow(irow);
    PlexPlaneId steelId = stRow -> GetPlaneId();

    std::string stName = GetGeoCompatibleName(steelId.AsString());
     
    if ( steelId.IsVetoShield() 
      && vldc.GetDetector() == Detector::kFar ) steelId =
      PlexVetoShieldHack::RenumberMuxToMdl(vldc,steelId);  // don't know why 
    else if ( steelId.GetPlane() > 485 ) continue; // because db has 0-497
    // Skip planes that have intentionally been placed at unphysical
    // locations
    if ( TMath::Abs(stRow->GetX0() ) > 50. ) continue;
    
    Double_t pairhalfx = 0;
    Double_t pairhalfy = 0;

    // If plane is first plane in supermodule, thickness of pair box is just
    // "Thickness", else pair box thickness is set by distance from back 
    // face of of current steel plane to back face of previous steel plane. 
    Int_t planeno = steelId.GetPlane();
    Double_t pairhalfz = 0;
    if ( vldc.GetDetector() == Detector::kCalDet || steelId.IsVetoShield() ) 
      pairhalfz = 0.5*(stRow->GetTotalZ())*fScale;
    else
      pairhalfz = 0.5*(stRow->GetThickness())*fScale;
    if ( !steelId.IsVetoShield() ) {
      if ( planeno != plnMin[0] ) {
        if ( nSM == 1 || (nSM > 1 && planeno != plnMin[1]) ) {
          const UgliDbiSteelPln* stPrevRow = steelTbl.GetRowByIndex(planeno-1);
          pairhalfz = 0.5*(stRow->GetZBack() - stPrevRow->GetZBack())*fScale;
        }
      }
    }
    
    if ( !steelId.IsVetoShield() ) {
      stplnVolume = new GeoSteelPlnVolume(this,steelId,
                                       stRow->GetThickness()*fScale);
      stplnVolume -> SetLineColor(kRed);
      stplnVolume -> SetVisibility(kTRUE);

      // Since we assume that pair is always placed at (x,y)=(0,0), need
      // to take into account bbox origin when determining pair dx,dy
      TGeoBBox* bbox = dynamic_cast<TGeoBBox*>(stplnVolume->GetShape());
      pairhalfx = bbox->GetDX() + TMath::Abs((bbox->GetOrigin())[0]);
      pairhalfy = bbox->GetDY() + TMath::Abs((bbox->GetOrigin())[1]);
    }
    
    // Scint plane
    GeoScintPlnVolume* scplnVolume = 0;
    const UgliDbiScintPln* scRow = 0;
    PlexPlaneId scintId = steelId;
    scintId.SetIsSteel(false);
    std::string scName = GetGeoCompatibleName(scintId.AsString());
    if ( PlaneCoverage::kNoActive != steelId.GetPlaneCoverage() && 
         PlaneCoverage::kUnknown  != steelId.GetPlaneCoverage() ) {
      scRow = ugliDbiTables.GetDbiScintPlnByIndex(scintId.GetPlane());
      MSG("Geo",Msg::kDebug) << "Build ScintPln " << scintId.AsString()
                             << " w/thickness " << scRow->GetThickness()*fScale
                             << endl;

      scplnVolume = new GeoScintPlnVolume(this,scintId,
                                          scRow->GetThickness()*fScale);
      scplnVolume -> SetVisibility(kFALSE); // deactivate when not debugging
      scplnVolume -> SetLineColor(kBlue);

      // This only works because scint planes are rotated in a way such that
      // width and height of bounding box don't change.
      TGeoBBox* bbox = dynamic_cast<TGeoBBox*>(scplnVolume->GetShape());
      pairhalfx = TMath::Max(pairhalfx,bbox->GetDX() 
                                     + TMath::Abs((bbox->GetOrigin())[0]));
      pairhalfy = TMath::Max(pairhalfy,bbox->GetDY() 
                                     + TMath::Abs((bbox->GetOrigin())[1]));
    }
    // avoid extension of near detector pair box into hall floor. Fix me.
    if ( vldc.GetDetector() == Detector::kNear ) pairhalfy = 2.91*fScale;
    
    // Build pair bounding box and position according to steel plane coords
    std::string pairName = GetGeoCompatibleName(steelId.AsString("b"));
    MSG("Geo",Msg::kDebug) << "Build PairBox " << pairName.c_str()
                           << " w/halfx,y,z " << pairhalfx << ","
                           << pairhalfy << "," << pairhalfz << "." << endl;

    Geo::EDetComponent detcomp = Geo::kPlnAirGap;
    if ( steelId.IsVetoShield() ) detcomp = Geo::kHall;
    TGeoVolume* pairVol = fGeoManager->MakeBox(pairName.c_str(),
                                             medMap.GetMedium(detcomp),
                                             pairhalfx,pairhalfy,pairhalfz);
    pairVol -> SetVisibility(kFALSE);
    pairVol -> VisibleDaughters(kTRUE);

    // Build 2 boxes representing air gaps.  "airgap0" is the gap between
    // the "front" (low-z face) of the scint plane in the current pair and the
    // "back" (high-z face) of the steel plane in the previous pair.
    // "airgap1" is the gap between the scint & steel plane in the current
    // pair.  The airgap coil segment will be added to each air gap.
    if ( !steelId.IsVetoShield() && coilVol ) {
      if ( planeno != plnMin[0] ) {
        if ( nSM == 1 || (nSM > 1 && planeno != plnMin[1]) ) {
          std::string airgap0 = pairName + "_air0";
          Double_t airgap0halfz = pairhalfz - 0.5*(stRow->GetTotalZ())*fScale;
          // totalz includes airgap even when scint pln not present? 
          // compensate for it here
          if ( !scplnVolume ) airgap0halfz 
                             = pairhalfz - 0.5*(stRow->GetThickness())*fScale;
          TGeoVolume* airgap0Vol = fGeoManager->MakeBox(airgap0.c_str(),
                                          medMap.GetMedium(Geo::kPlnAirGap),
                                          pairhalfx,pairhalfy,airgap0halfz);
          airgap0Vol -> SetVisibility(kFALSE);
          airgap0Vol -> VisibleDaughters(kTRUE);
          Double_t airgap0z0 = -pairhalfz + airgap0halfz;
          pairVol->AddNode(airgap0Vol,1,new TGeoTranslation(0.,0.,airgap0z0));
          if (dettype == Detector::kNear) airgap0Vol->AddNode(coilVol,1,rot45);
          else airgap0Vol -> AddNode(coilVol,1,gGeoIdentity);
        }
      }
      if ( scplnVolume ) {  // can't have airgap1 without active plane
        std::string airgap1 = pairName + "_air1";
        Double_t airgap1halfz = 0.5*(stRow->GetTotalZ()-stRow->GetThickness() 
                          - scRow->GetThickness())*fScale;
        TGeoVolume* airgap1Vol = fGeoManager->MakeBox(airgap1.c_str(),
                                     medMap.GetMedium(Geo::kPlnAirGap),
                                     pairhalfx,pairhalfy,airgap1halfz);
        airgap1Vol -> SetVisibility(kFALSE);
        airgap1Vol -> VisibleDaughters(kTRUE);
        Double_t airgap1z0 = pairhalfz-stRow->GetThickness()*fScale
                           - airgap1halfz;
        pairVol -> AddNode(airgap1Vol,1,new TGeoTranslation(0.,0.,airgap1z0));

        if ( dettype == Detector::kNear ) airgap1Vol->AddNode(coilVol,1,rot45);
        else airgap1Vol -> AddNode(coilVol,1,gGeoIdentity);
      }  
    }
    
    // Add to shield groups so as to be able to build bounding volume
    // default for caldet shield planes is still to place pair directly in hall
    std::string pairPath = fHallPath + "/" + pairName +"_1"; 
    if ( steelId.IsVetoShield() && vldc.GetDetector() == Detector::kFar ) { 
      GeoShield::EGroupType group=fGeoShield.AddVolume(pairVol,steelId,stRow);
      pairPath = fHallPath + "/" + GeoShield::AsString(group) + "_1/" + 
                 pairName+"_1";
    }
      
    int iSM = -1; // veto shield
    if ( !steelId.IsVetoShield() ) {
      iSM = 0;
      if ( nSM > 1 && steelId.GetPlane() > plnMax[0] ) iSM = 1;
    }
      
    if ( stplnVolume ) {
      // Place the steel plane within the pairVol
      std::string stplnNodeName = stName + "_1";
      std::string steelPlnPath = pairPath + "/" + stplnNodeName;
      double z0steel = pairhalfz - 0.5*(stRow->GetThickness())*fScale;
      MSG("Geo",Msg::kDebug) << "Placing steel pln at " << 0. << ","
                             << 0. << "," << z0steel << endl;
      GeoSteelPlnNode* stplnNode = new GeoSteelPlnNode(this,stplnVolume,
                          new TGeoTranslation(stName.c_str(),0.,0.,z0steel),
                          pairVol,steelPlnPath,stplnNodeName,steelId);
      
      // fSMZMin/Max in MARS coordinates
      if ( vldc.GetDetector() != Detector::kCalDet ) 
        fSMZMin[iSM] = TMath::Min(fSMZMin[iSM],(stRow->GetZBack()
                                  -stRow->GetThickness())*fScale);  
      else 
        fSMZMin[iSM] = TMath::Min(fSMZMin[iSM],(stRow->GetZBack()
                                  -stRow->GetTotalZ())*fScale);      
      fSMZMax[iSM] = TMath::Max(fSMZMax[iSM],(stRow->GetZBack())*fScale);

      // Insert steelplnNode in plane map
      fPlaneMap[steelId] = stplnNode;
    }
    
    if ( scplnVolume ) {
      // Place the scint plane within the pairVol
      Float_t x0pln = (scRow->GetX0RelSteel())*fScale;
      Float_t y0pln = (scRow->GetY0RelSteel())*fScale;
      Float_t z0pln = pairhalfz + 
                      (-stRow->GetTotalZ()+0.5*(scRow->GetThickness()))*fScale;

      // Rotation matrix relative to steel
      // (thetax,phix) = (90, 0+zrotrelsteel)
      // (thetay,phiy) = (90,90+zrotrelsteel)
      // (thetaz,phiz) = ( 0, 0)
      Double_t zrotrelsteel = scRow->GetZRotRelSteelDeg();

      MSG("Geo",Msg::kDebug) << "Placing scint pln at " << x0pln << ", "
                             << y0pln << "," << z0pln 
                             << " zrotrelsteel " << zrotrelsteel << endl;
        
      TGeoRotation* scRotMatrix = new TGeoRotation(scName.c_str(),90.,
                                  zrotrelsteel,90.,90.+zrotrelsteel,0.,0.);
      scRotMatrix -> RegisterYourself();
      
      TGeoCombiTrans* combiMatrix 
         = new TGeoCombiTrans(scName.c_str(),x0pln,y0pln,z0pln,scRotMatrix);

      std::string scplnNodeName = scName + "_1";
      std::string scintPlnPath = pairPath + "/" + scplnNodeName;
      GeoScintPlnNode* scplnNode 
          = new GeoScintPlnNode(this,scplnVolume,combiMatrix,pairVol,
                                scintPlnPath,scplnNodeName,scintId);
    
      this -> BuildModules(ugliDbiTables, scplnNode);

      // Insert scint pln node in plane map
      fPlaneMap[scintId] = scplnNode;

    }
    
  } // end of loop over all plane db rows

  
  return;
  
}
  
//_____________________________________________________________________________
Int_t GeoGeometry::GetSMPlaneLimits(const VldContext& vldc,
                                    Int_t* plnMin, Int_t* plnMax) {
  // Private, but could be public?
  // Determine Super Module plane limits for specified vldc.  
  // Returns number of super modules.  plnMin/Max[0] is filled for
  // supermodule 0 and plnMin/Max[1] for supermodule 1 if applicable.

  Detector::Detector_t dettype = vldc.GetDetector();
  Int_t nSM = 0;
  switch (dettype) {
  case Detector::kFar: 
    nSM = 2;
    plnMin[0] = 0;
    plnMax[0] = 248;
    plnMin[1] = 249;
    plnMax[1] = 485;
    break;
  case Detector::kNear:
    nSM = 1;
    plnMin[0] = 0;
    plnMax[0] = 281;
    break;
  case Detector::kCalDet:
    nSM = 1;
    plnMin[0] = 0;
    plnMax[0] = 60; // main detector block, 61-64 on floor
    break;
  default:
    MSG("Geo",Msg::kError) 
    << "GetSMPlaneLimits failed for unknown detector type\n"
    << Detector::AsString(dettype) << "." << endl;
  } // end of switch

  return nSM;
  
}

//_____________________________________________________________________________
TGeoVolume* GeoGeometry::BuildFarMARS() {
  // Private. Build MARS appropriate for the far detector site.

  UpdateGlobalManager();
 
  MSG("Geo",Msg::kDebug) << "GeoGeometry::BuildFarMARS " << endl;
  assert(fGeoManager);

  const GeoMediumMap& medMap = GetMediumMap();

  // Dimensions of MARS should be db parameters?
  Float_t mars_halfx =  85.*fScale; //  85 m MARS half-x
  Float_t mars_halfy =  30.*fScale; //  30 m MARS half-y
  Float_t mars_halfz = 250.*fScale; // 250 m MARS half-z
  
  // Far detector rock medium is different than that of near
  TGeoVolume* volMARS = fGeoManager->MakeBox("MARS",
                                        medMap.GetMedium(Geo::kMars),
                                        mars_halfx,mars_halfy,mars_halfz);

  // Build shaft, tunnel, Soudan 2 cavern & ghost of S2 detector out of air 
  // and according to gminos dimensions, except as noted.

  Float_t shft_halfx =  1.*fScale; //  1 m shaft half-x
  Float_t shft_halfy = 30.*fScale; // 30 m shaft half-y == MARS half-y
  Float_t shft_halfz =  1.*fScale; //  1 m shaft half-z
  TGeoVolume* volSHFT = gGeoManager->MakeBox("SHFT",
                                             medMap.GetMedium(Geo::kHall),
                                             shft_halfx,shft_halfy,shft_halfz);

  // Shorten tunnel by 1 m from gminos dimensions to avoid overlap with LINR
  Float_t tunl_halfx =  2.*fScale; //  2 m tunnel half-x
  Float_t tunl_halfy =  2.*fScale; //  2 m tunnel half-y 
  Float_t tunl_halfz = 13.5*fScale; // 13.5 m tunnel half-z
  TGeoVolume* volTUNL = gGeoManager->MakeBox("TUNL",
                                            medMap.GetMedium(Geo::kHall),
                                            tunl_halfx,tunl_halfy,tunl_halfz);

  Float_t s2hl_halfx =  7.  *fScale; //  7    m S2 hall half-x
  Float_t s2hl_halfy =  5.65*fScale; //  5.65 m S2 hall half-y 
  Float_t s2hl_halfz = 35.  *fScale; // 35.   m S2 hall half-z
  TGeoVolume* volS2HL = gGeoManager->MakeBox("S2HL",
                                            medMap.GetMedium(Geo::kHall),
                                            s2hl_halfx,s2hl_halfy,s2hl_halfz);

  // The ghost of the S2 detector size is adjusted from that in gminos to 
  // be more realistic, but it doesn't matter of course because it's not there!
  Float_t s2dt_halfx =  4. *fScale; // 4   m S2 detector half-x
  Float_t s2dt_halfy =  2.7*fScale; // 2.7 m S2 detector half-y (4.5, gminos) 
  Float_t s2dt_halfz =  7.5*fScale; // 7.5 m S2 detector half-z 
  TGeoVolume* volS2DT = gGeoManager->MakeBox("S2DT",
                                            medMap.GetMedium(Geo::kHall),
                                            s2dt_halfx,s2dt_halfy,s2dt_halfz);
  volS2HL -> AddNode(volS2DT,1,
                     new TGeoTranslation(0*fScale,-1.15*fScale,-23.5*fScale));


  volMARS -> AddNode(volSHFT,1,
                     new TGeoTranslation(-5.*fScale,0.*fScale,98.18*fScale));
  volMARS -> AddNode(volTUNL,1,
                     new TGeoTranslation(-5.*fScale,-4.1*fScale,83.68*fScale));
  TGeoRotation* s2rot = new TGeoRotation("s2rot",64,0,90,90,26,180);
  s2rot -> RegisterYourself();
  volMARS -> AddNode(volS2HL,1,
                     new TGeoCombiTrans(36.*fScale,-0.45*fScale,61.68*fScale,
                                        s2rot));

  return volMARS;

}

//_____________________________________________________________________________
TGeoVolume* GeoGeometry::BuildNearMARS() {
  // Private. Build MARS appropriate for the near detector site.

  UpdateGlobalManager();
 
  MSG("Geo",Msg::kDebug) << "GeoGeometry::BuildNearMARS " << endl;
  assert(fGeoManager);

  const GeoMediumMap& medMap = GetMediumMap();

  // Dimensions of MARS should be db parameters?
  Float_t mars_halfx =  85.*fScale; //  85 m MARS half-x
  Float_t mars_halfy =  30.*fScale; //  30 m MARS half-y
  Float_t mars_halfz = 250.*fScale; // 250 m MARS half-z
  
  // Near detector rock medium is different than that of far 
  TGeoVolume* volMARS = fGeoManager->MakeBox("MARS",
                                           medMap.GetMedium(Geo::kMars),
                                           mars_halfx,mars_halfy,mars_halfz);

  return volMARS;

}

//_____________________________________________________________________________
void GeoGeometry::BuildNearCoil(TGeoVolume* hallVol) {
  // Private. Build the near detector coil. Fix me.  
  // x0,y0 is coil hole center in hall coordinates (assumed constant
  // over all planes).

  UpdateGlobalManager();
 
  MSG("Geo",Msg::kDebug) << "GeoGeometry::BuildNearCoil " << endl;
  assert(fGeoManager);

  const GeoMediumMap& medMap = GetMediumMap();

  // The dimensions of the coil should be extracted from the db
  Float_t shalfx  = Geo::kNeckRad*fScale; // steel 
  Float_t shalfy  = shalfx; // steel
  Float_t ahalfx  = Geo::kThroatRad*fScale; // air
  Float_t ahalfy  = ahalfx; // air
  Float_t alhalfx = Geo::kCoilRad*fScale; // aluminum
  Float_t alhalfy = alhalfx; // aluminum
 
  // Build coil along z-direction, extends to ends of super module
  Float_t zmin,zmax;
  GetZExtent(zmin,zmax,-1); // returned in MARS coordinates
  // Convert to HALL coordinates
  TGeoNode* linrNode = fGeoManager->GetVolume("MARS")->GetNode("LINR_1");
  const Double_t* h0 = linrNode -> GetMatrix() -> GetTranslation();
  zmin -= h0[2];
  zmax -= h0[2];
  Float_t x0 = -h0[0];
  Float_t y0 = -h0[1];
  
  Float_t zcoillen = zmax - zmin;
  TGeoVolume* BUS0 = fGeoManager->MakeBox("BUS0",
                                           medMap.GetMedium(Geo::kRCShaft),
                                           shalfx,shalfy,0.5*zcoillen);
  BUS0 -> SetLineColor(14);  // gray  
  BUS0 -> SetVisibility(kTRUE);

  TGeoVolume* BUA0 = fGeoManager->MakeBox("BUA0",
                                           medMap.GetMedium(Geo::kRCAir),
                                           ahalfx,ahalfy,0.5*zcoillen);
  BUA0 -> SetLineColor(14);  
  BUA0 -> SetVisibility(kTRUE);
  BUS0 -> AddNode(BUA0,1,gGeoIdentity);

  TGeoVolume* BUC0 = fGeoManager->MakeBox("BUC0",
                                           medMap.GetMedium(Geo::kRCCoil),
                                           alhalfx,alhalfy,0.5*zcoillen);
  BUC0 -> SetLineColor(14);  
  BUC0 -> SetVisibility(kTRUE);
  BUA0 -> AddNode(BUC0,1,gGeoIdentity);

  // Rotate first, then translate
  TGeoRotation *rot45 = GetRot45();
  Float_t zcoil = (zmin+zmax)/2.;

  // Place return coil in hall
  Float_t steelhalfy = 1.9*fScale; // 1.9 m halfwidth in y direction
  Float_t xcoil = x0 - steelhalfy;
  Float_t ycoil = y0 - steelhalfy;
  // overlap because it overlaps pair plane volumes. Fix me.
  hallVol -> AddNodeOverlap(BUS0,2,
                            new TGeoCombiTrans(xcoil,ycoil,zcoil,rot45)); 

  // Along xy face
  Float_t cos45 = 0.707107;
  Float_t coillenxy = steelhalfy/cos45 + 2.*shalfx;
  TGeoVolume* BUSXY = fGeoManager->MakeBox("BUSXY",
                                            medMap.GetMedium(Geo::kRCShaft),
                                            shalfx,0.5*coillenxy,shalfy);
  BUSXY -> SetLineColor(14);  // gray  
  BUSXY -> SetVisibility(kTRUE);
  TGeoVolume* BUAXY = fGeoManager->MakeBox("BUAXY",
                                            medMap.GetMedium(Geo::kRCAir),
                                            ahalfx,0.5*coillenxy,ahalfy);
  BUAXY -> SetLineColor(14);  
  BUAXY -> SetVisibility(kTRUE);
  BUSXY -> AddNode(BUAXY,1,gGeoIdentity);

  TGeoVolume* BUCXY = fGeoManager->MakeBox("BUCXY",
                                            medMap.GetMedium(Geo::kRCCoil),
                                            alhalfx,0.5*coillenxy,alhalfy);
  BUCXY -> SetLineColor(14);  
  BUCXY -> SetVisibility(kTRUE);
  BUAXY -> AddNode(BUCXY,1,gGeoIdentity);

  // Will be rotated and then translated
  TGeoRotation* rotneg45 = GetRot315();
  xcoil = x0 - (0.5*coillenxy - shalfx)*cos45;
  ycoil = y0 - (0.5*coillenxy - shalfy)*cos45;
  zcoil = zmin - shalfy;
  hallVol -> AddNode(BUSXY,1,new TGeoCombiTrans(xcoil,ycoil,zcoil,rotneg45));  

  zcoil = zmax + shalfy;
  hallVol -> AddNode(BUSXY,2,new TGeoCombiTrans(xcoil,ycoil,zcoil,rotneg45));  

  return;

}

//_____________________________________________________________________________
void GeoGeometry::BuildFarCoil(TGeoVolume* hallVol) {
  // Private. Build the far detector coil outside of the detector.  
  // x0,y0 is coil hole center in HALL coordinates (assumed constant
  // over all planes).

  UpdateGlobalManager();
 
  MSG("Geo",Msg::kDebug) << "GeoGeometry::BuildFarCoil " << endl;
  assert(fGeoManager);

  const GeoMediumMap& medMap = GetMediumMap();

  // These dimensions/etc. of the coil should be in the db, but for now
  // use hardwired values

  // Nomenclature (BUS) is chosen to match that of gminos
  // BUSn is the center of the coil made up of medium FARCOIL
  // In the convention used here, n = supermodule number

  Float_t sRad = Geo::kCoilRad*fScale; // coil radius 

  Float_t halfhally = dynamic_cast<TGeoBBox*>(hallVol -> GetShape())->GetDY();
  
  // SM0 coil along z-direction of detector, ends at supermodule ends
  // return coil lies on floor
  Float_t zmin0,zmax0;
  GetZExtent(zmin0,zmax0,0);
  // Convert to HALL coordinates
  TGeoNode* linrNode = fGeoManager->GetVolume("MARS")->GetNode("LINR_1");
  const Double_t* h0 = linrNode -> GetMatrix() -> GetTranslation();
  zmin0 -= h0[2];
  zmax0 -= h0[2];
  Float_t x0 = -h0[0];
  Float_t y0 = -h0[1];

  // Per specification from J. Nelson received 8/1/2006, model return
  // coil lying along far detector floor as box of halfx 0.225 m and
  // halfy 0.05 m.
  Float_t retcoilhalfx = 0.225*fScale;
  Float_t retcoilhalfy = 0.05*fScale;

  Float_t zcoillen0 = zmax0-zmin0;
  TGeoVolume* BUS0=fGeoManager->MakeBox("BUS0",medMap.GetMedium(Geo::kRCCoil),
                                      retcoilhalfx,retcoilhalfy,0.5*zcoillen0);
  BUS0 -> SetLineColor(42); // tan  
  BUS0 -> SetVisibility(kTRUE);

  // Place SM0 return coil in hall
  Float_t z0 = 0.5*(zmin0 + zmax0);
  hallVol -> AddNode(BUS0,1,new TGeoTranslation("BUS0",x0,
                                           -halfhally+retcoilhalfy,z0));
  
  // SM1 coil along z-direction of detector, ends at supermodule ends
  Float_t zmin1,zmax1;
  GetZExtent(zmin1,zmax1,1);
  // Convert to HALL coordinates
  zmin1 -= h0[2];
  zmax1 -= h0[2];

  Float_t zcoillen1 = zmax1-zmin1;
  TGeoVolume* BUS1=fGeoManager->MakeBox("BUS1",medMap.GetMedium(Geo::kRCCoil),
                                      retcoilhalfx,retcoilhalfy,0.5*zcoillen1);
  BUS1 -> SetLineColor(42);  // tan  
  BUS1 -> SetVisibility(kTRUE);

  // Place SM1 return coil in hall
  Float_t z1 = 0.5*(zmin1 + zmax1);
  hallVol -> AddNode(BUS1,1,new TGeoTranslation("BUS1",x0,
                                                -halfhally+retcoilhalfy,z1));

  // Build coil volume on xy faces
  // Dimensions for horse's tail are roughly those received from J. Nelson 
  // correspondence of 8/1/2006, but y-length of bottom section is adjusted to 
  // actual distance to floor, and top & mid-sections are trapezoids 
  // instead of boxes, dimension in z has been squeezed to fit the supermodule
  // gap...
  Float_t ycoillen = y0+halfhally+sRad; // top of coil to floor

  // gap between bottom section of tail & plane. Adjust to fit avail region.
  Float_t zbotgap = 0.10*fScale; 
  // Define bounding box to contain tail 
  Float_t xhalftail = 0.40*fScale;
  Float_t yhalftail = 0.5*ycoillen;
  Float_t zhalftail = 0.5*(zbotgap+0.4*fScale);
  TGeoVolume* BUXY = fGeoManager->MakeBox("BUXY",medMap.GetMedium(Geo::kRCAir),
                                           xhalftail,yhalftail,zhalftail);
  BUXY -> SetLineColor(14); // gray bounding box, invisible by default 
  BUXY -> SetVisibility(kFALSE);

  // Horse's tail bounding box is placed at 4 end planes, rotating &
  // translating as necessary
  // rotFlip defines a 180 rotation
  TGeoRotation *rotFlip = new TGeoRotation("BUXY",90,180,90,90,180,0);
  rotFlip -> RegisterYourself();
  
  Float_t xcoil = x0;
  Float_t ycoil = y0 + sRad - ycoillen/2.;
  hallVol -> AddNode(BUXY,1,new TGeoTranslation("BUXY_1",xcoil,ycoil,
                     zmin0-zhalftail));
  hallVol -> AddNode(BUXY,2,new TGeoCombiTrans("BUXY_2",xcoil,ycoil,
                     zmax0+zhalftail,rotFlip));
  hallVol -> AddNode(BUXY,3,new TGeoTranslation("BUXY_3",xcoil,ycoil,
                     zmin1-zhalftail));
  hallVol -> AddNode(BUXY,4,new TGeoCombiTrans("BUXY_4",xcoil,ycoil,
                     zmax1+zhalftail,rotFlip));

  // Now build the inserts for the coil tail bounding box
  // Top section is a trapezoid
  Float_t yhalftop = 0.5*0.53*fScale;
  Float_t zhalftop = 0.10*fScale;
  Float_t xhalfloy = 0.2825*fScale;
  Float_t xhalfhiy = 0.12*fScale;
  TGeoTrap* bxyTopShp = new TGeoTrap("BXYTOP",zhalftop,0.,0.,
                                     yhalftop,xhalfloy,xhalfhiy,0.,
                                     yhalftop,xhalfloy,xhalfhiy,0.);
  TGeoVolume* BXYTOP = new TGeoVolume("BXYTOP",bxyTopShp,
                                       medMap.GetMedium(Geo::kRCCoil));
  BXYTOP -> SetLineColor(42);  // tan
  BXYTOP -> SetVisibility(kTRUE);
  BUXY -> AddNode(BXYTOP,1,new TGeoTranslation(0.,yhalftail-sRad-yhalftop,
                                              -zhalftail+zhalftop));
  

  // Middle section is a trapezoid that skews towards the face of the
  // steel plane with decreasing y
  Float_t yhalfmid = 0.5*0.85*fScale;
  TGeoArb8* bxyMidShp = new TGeoArb8(yhalfmid);
  bxyMidShp -> SetVertex(0,-0.4*fScale,0.05*fScale);
  bxyMidShp -> SetVertex(1,-0.4*fScale,0.15*fScale);
  bxyMidShp -> SetVertex(2,0.4*fScale,0.15*fScale);
  bxyMidShp -> SetVertex(3,0.4*fScale,0.05*fScale);
  bxyMidShp -> SetVertex(4,-0.2825*fScale,-0.25*fScale);
  bxyMidShp -> SetVertex(5,-0.2825*fScale,-0.05*fScale);
  bxyMidShp -> SetVertex(6,0.2825*fScale,-0.05*fScale);
  bxyMidShp -> SetVertex(7,0.2825*fScale,-0.25*fScale);
  
  TGeoVolume* BXYMID = new TGeoVolume("BXYMID",bxyMidShp,
                                       medMap.GetMedium(Geo::kRCCoil));
  BXYMID -> SetLineColor(42); // tan
  BXYMID -> SetVisibility(kTRUE);

  TGeoRotation *rotMid = new TGeoRotation("BXYMid",90,0,0,0,90,90);
  rotMid -> RegisterYourself();
  BUXY -> AddNode(BXYMID,1,new TGeoCombiTrans("BXYMID",0.,
                  yhalftail-sRad-2.*yhalftop-yhalfmid,0.,rotMid));

  // Bottom section is a box that extends to floor
  Float_t xhalfbot = 0.4*fScale;
  Float_t yhalfbot = 0.5*(ycoillen-sRad-2.*yhalftop-2.*yhalfmid);
  Float_t zhalfbot = 0.05*fScale;
  
  TGeoVolume* BXYBOT = fGeoManager -> MakeBox("BXYBOT",
                                               medMap.GetMedium(Geo::kRCCoil),
                                               xhalfbot,yhalfbot,zhalfbot);
  BXYBOT -> SetLineColor(42); // tan
  BXYBOT -> SetVisibility(kTRUE);
  BUXY -> AddNode(BXYBOT,1,new TGeoTranslation(0.,-yhalftail+yhalfbot,
                                               zhalftail-zbotgap-zhalfbot));
  
  // Return coil along floor that extends outside of detector
  TGeoVolume* BXYRET = fGeoManager -> MakeBox("BXYRET",
                               medMap.GetMedium(Geo::kRCCoil),
                               retcoilhalfx,retcoilhalfy,0.5*zbotgap);
  BXYRET -> SetLineColor(42);  // tan
  BXYRET -> SetVisibility(kTRUE);
  BUXY -> AddNode(BXYRET,1,new TGeoTranslation(0.,-yhalftail+retcoilhalfy,
                                               zhalftail-0.5*zbotgap));

  // Segment of coil tube that extends out beyond detector
  TGeoVolume* BXYSTB = fGeoManager -> MakeTube("BXYSTB",
                                                medMap.GetMedium(Geo::kRCCoil),
                                                0.,sRad,zhalftail);
  BXYSTB -> SetLineColor(42);  // tan
  BXYSTB -> SetVisibility(kTRUE);
  BUXY -> AddNodeOverlap(BXYSTB,1,new TGeoTranslation(0.,yhalftail-sRad,0.));
                                   
  return;

}

//_____________________________________________________________________________
void GeoGeometry::BuildModules(const UgliDbiTables& ugliDbiTables,
                               GeoScintPlnNode* plnNode) {
  // Private. Purpose: Build the modules as part of the scintillator plane
  
  UpdateGlobalManager();
  
  assert(fGeoManager);

  const PlexPlaneId& plexPlaneId = plnNode -> GetPlexPlaneId();
  MSG("Geo",Msg::kDebug) <<"GeoGeometry::BuildModules for scint pln  " 
                         << plnNode->GetName() << endl;

  // Global path to plane node. Used to build module global path.
  std::string plnPath = plnNode->GetGlobalPath();

  // ScintPlnStruct used to retrieve number of modules in plane  
  UgliDbiStructHash planestructhash(plexPlaneId);
  const UgliDbiScintPlnStruct* plnStruct = ugliDbiTables.fScintPlnStructTbl
        .GetRowByIndex(planestructhash.HashAsPlane());
  Int_t nmodules = plnStruct -> GetNModules();
  
  // Loop over all modules in plane
  for ( int imod = 0; imod < nmodules; imod++ ) {
    // For each module
    PlexScintMdlId plexModuleId(plexPlaneId,imod);

    std::string mdlName = GetGeoCompatibleName(plexModuleId.AsString());
    
    // Build the module volume
    // Ideally could just build one module of each type but alignment varies
    // tpos and this may cause slight size variations
    // So for now generate one mdl volume per module
    GeoScintMdlVolume* mdlVol = new GeoScintMdlVolume(this,plexModuleId,
                                                      ugliDbiTables);
 
    // Build the node and place it in the plane
    // Establish the rotational & translational matrix for the module
    const UgliDbiScintMdl* scModule 
      = ugliDbiTables.GetDbiScintMdlById(plexModuleId);
    Float_t tposMod = (scModule -> GetTPosRelPln())*fScale;
    Float_t lposMod = (scModule -> GetLPosRelPln())*fScale;
    Float_t zrotMod = scModule -> GetZRotRelPlnDeg(); 
 
    TGeoRotation* modRotMatrix = new TGeoRotation(mdlName.c_str(), 90.,zrotMod,
                                                  90.,90.+zrotMod,0.,0.);
    modRotMatrix -> RegisterYourself();

    TGeoCombiTrans* combiMatrix = new TGeoCombiTrans(mdlName.c_str(),
                                     lposMod,tposMod,0.,modRotMatrix);
    std::string mdlNodeName = mdlName + "_1";

    MSG("Geo",Msg::kVerbose) << "Mdl node " << mdlNodeName.c_str() 
                             << "tpos, lpos, zrot(deg) rel pln " 
                             << tposMod << "," << lposMod << "," << zrotMod 
                             << endl;
    
    GeoScintMdlNode* mdlNode = new GeoScintMdlNode(this,mdlVol,combiMatrix,
                                                   plnNode,mdlNodeName,
                                                   plexModuleId,
                                          scModule->GetClearLenEast()*fScale,
                                          scModule->GetClearLenWest()*fScale,
                                          scModule->GetWlsLenEast()*fScale,
                                          scModule->GetWlsLenWest()*fScale);

    // Add strips to module node
    BuildStrips(ugliDbiTables,mdlNode);

  }
  
}
  
//_____________________________________________________________________________
void GeoGeometry::BuildStrips(const UgliDbiTables& ugliDbiTables,
                              GeoScintMdlNode* mdlNode) {
  // Purpose:: Build this planes's scint strips and add them as nodes to the
  // det plane volume
  // This method should be moved to GeoScintPlnVolume

  UpdateGlobalManager();
  assert(fGeoManager);
 
  PlexScintMdlId plexModuleId = mdlNode->GetPlexScintMdlId();
  MSG("Geo",Msg::kDebug) << "GeoGeometry::BuildStrips for module " 
                         << mdlNode->GetName() << endl;
  
  // Determine lo,hi strip in plane
  UgliDbiStructHash modulestructhash(plexModuleId);
  const UgliDbiScintMdlStruct* scModuleStruct 
       = ugliDbiTables.fScintMdlStructTbl
        .GetRowByIndex(modulestructhash.HashAsScintMdl());
  Int_t lostrip = scModuleStruct -> GetFirstStrip();
  Int_t histrip = scModuleStruct -> GetLastStrip();

  for ( int istp = lostrip; istp <= histrip; istp++ ) {
    // For each strip in module
    PlexStripEndId seid(plexModuleId,istp);
    std::string stpName = GetGeoCompatibleName(seid.AsString());

    MSG("Geo",Msg::kVerbose)  << "Strip " << istp << " name " 
                              << stpName.c_str() << endl;
    const UgliDbiStrip* scStrip = ugliDbiTables.GetDbiStripById(seid);

    // The strip volumes have been pre-built
    std::string stpVolName = GetGeoCompatibleName(seid.AsString("s"));
    TGeoVolume* stpVol = fGeoManager -> GetVolume(stpVolName.c_str());
    if ( !stpVol ) {
      MSG("Geo",Msg::kFatal) 
         << "Unable to find pre-built stripvol for volume "
         << stpName.c_str() << " shape name " 
         << stpVolName.c_str() << endl;
      abort();
    } 

    // Build the node and place it in the plane
    Float_t tposStp = (scStrip -> GetTPosRelMdl())*fScale;
    Float_t lposStp = (scStrip -> GetLPosRelMdl())*fScale;
    Float_t zrotStp = scStrip -> GetZRotRelMdlDeg();
    MSG("Geo",Msg::kVerbose) << "tpos, lpos, zrot(deg) rel to module " 
                             << tposStp << ", " << lposStp << ", " 
                             << zrotStp << endl;

    std::string stpNodeName = stpVolName+"_1";

    TGeoRotation* stpRotMatrix = new TGeoRotation(stpNodeName.c_str(), 90.,
                                       zrotStp,90.,90.+zrotStp,0.,0.);
    stpRotMatrix -> RegisterYourself();
    TGeoCombiTrans* combiMatrix 
     = new TGeoCombiTrans(stpNodeName.c_str(),lposStp,tposStp,0.,stpRotMatrix);

    new GeoStripNode(this,stpVol,combiMatrix,mdlNode,stpNodeName,seid);

  }

  return;
  
}

//__________________________________________________________________________ 
void GeoGeometry::BuildStripVolumes(const VldContext& vldc) {
  // pre-build set of strip volumes for this vldc.  strip volumes will be 
  // placed appropriately in module volumes as nodes when planes are 
  // constructed.

  UpdateGlobalManager();

  MSG("Geo",Msg::kDebug) << "GeoGeometry::BuildStripVolumes " << endl;
  assert(fGeoManager);

  DbiResultPtr<UgliDbiStripStruct> stripTbl(vldc);
  Detector::Detector_t dettype = vldc.GetDetector();
  
  for ( unsigned int irow = 0; irow < stripTbl.GetNumRows(); irow++ ) {
    const UgliDbiStripStruct* sRow = stripTbl.GetRow(irow);
    string stpName = GetGeoCompatibleName(sRow -> GetShapeName());
    PlaneView::PlaneView_t plnview = sRow->GetPlaneView();
    bool isVetoShield = false;
    if ( dettype == Detector::kFar ) {
      if ( plnview != PlaneView::kU && plnview != PlaneView::kV ) 
        isVetoShield = true;
    }
    else if ( dettype == Detector::kCalDet ) {
      if ( plnview != PlaneView::kX && plnview != PlaneView::kY )
        isVetoShield = true;
    }
    
    // Volumes will self-register with fGeoManager
    Float_t totallen = (sRow->GetTotalLen())*fScale;
    Float_t leneastpart = (sRow->GetLenEastPart())*fScale;
    Float_t lenwestpart = (sRow->GetLenWestPart())*fScale;
    // Clean up some caldet database crud
    if ( totallen < 1.E-4*fScale ) continue;
    if ( leneastpart < 1.E-4*fScale ) leneastpart = totallen;
    if ( lenwestpart < 1.E-4*fScale ) lenwestpart = totallen;
    
    new GeoStripVolume(this,stpName.c_str(),totallen,leneastpart,lenwestpart,
                                       (sRow->GetWlsLenEast())*fScale, 
                                       (sRow->GetWlsLenWest())*fScale, 
                                       (sRow->GetWlsLenBypass())*fScale,
                                       isVetoShield);
  } 

}

//_____________________________________________________________________________
void GeoGeometry::UpdateNodeMatrices(TGeoVolume* volume) {
  // Private method. Loop through all GeoNodes in geometry starting at volume 
  // and construct global matrices for use in local to global (MARS) 
  // coordinates

  UpdateGlobalManager();
  if (!volume) {
    MSG("Geo",Msg::kFatal) << "UpdateNodeMatrices called with null volume ptr "
                           << endl;
    abort();
  }
  
  Int_t ndaughter = volume -> GetNdaughters();
  for ( int id = 0; id < ndaughter; id++ ) {
    TGeoNode* daughternode = volume->GetNode(id);
    GeoNode* geonode = dynamic_cast<GeoNode*>(daughternode);
    if ( geonode ) geonode -> UpdateGlobalMatrix();
    
    TGeoVolume* daughtervol = daughternode->GetVolume();
    this -> UpdateNodeMatrices(daughtervol);
  }
  
}

//_____________________________________________________________________________
TGeoVolume* GeoGeometry::BuildCoilAirGapVolume(const VldContext& vldc) const {
  // Build coil region volume to be used in airgaps between scint and
  // steel planes.

  UpdateGlobalManager();
  
  TGeoVolume* coilVol = 0;

  Detector::Detector_t dettype = (Detector::Detector_t)vldc.GetDetector();

  switch (dettype) {

  case Detector::kFar:
    coilVol = BuildFarCoilAirGapVolume();
    break;

  case Detector::kNear:
    coilVol = BuildNearCoilAirGapVolume();
    break;

  default:
    break;

  }

  return coilVol;
  
}


//_____________________________________________________________________________
TGeoVolume* GeoGeometry::BuildFarCoilAirGapVolume() const {
  //
  // Build far coil region to be used in airgaps between scint and steel
  // planes.  The far detector coil in the air gap
  // regions is built according to the same parameters as are used in
  // the build of the coil in the far detector scint pln:
  //   1)Insert a "bypass" (r=Geo::kFarBypassRad) region (not to be confused
  //     with the bypass plastic insert).  This is constructed of magnetized
  //     air, and indicates the region near the coil where the detailed
  //     field map will be used.
  //   2)Insert into the "bypass" a "flange" (r=Geo::kFlangeRad) of magnetized
  //     iron.
  //   3)Insert into the "flange" a "throat" (r=Geo::kThroatRad) of magnetized
  //     air.
  //   4)Insert into the "throat" a "coil segment" (r=Geo::kCoilRad) of
  //     magnetized FARCOIL medium.  The coil segment is set to rest on the
  //     throat (off-center in y).
  //
  // The thickness in z of all pieces is set to be -1 when defining the shapes,
  // and built as a TGeoVolumeMulti.  When placed as a node in the air gap
  // volume, the volumemulti shape will take on the thickness of its parent.
  // This allows for variable thickness of the air gaps in a non-perfect
  // geometry.
  //

  UpdateGlobalManager();

  // The nomenclature used here is "F" (far), "G" (gap), followed by
  // the two letter code to indicate the coil part type.

  const GeoMediumMap& medMap = GetMediumMap();
  
  // First check to see if volumemulti has been built. Need only build once.
  TGeoVolume* volBypass = gGeoManager->GetVolume("FGBY");
  if ( volBypass ) return volBypass;

  Double_t scale = Geo::GetScale(GetAppType());

  TGeoVolume* volCoil
    = gGeoManager->MakeTube("FGCO",medMap.GetMedium(Geo::kCRGapCoil),
                            0,Geo::kCoilRad*scale,-1.); // coil
  volCoil -> SetLineColor(42); // tan

  TGeoVolume* volThroat
    = gGeoManager->MakeTube("FGTR",medMap.GetMedium(Geo::kCRGapThroat),
                            0,Geo::kThroatRad*scale,-1.); // throat
  volThroat -> SetLineColor(kBlue);

  TGeoVolume* volFlange
    = gGeoManager->MakeTube("FGFL",medMap.GetMedium(Geo::kCRGapFlange),
                             0,Geo::kFlangeRad*scale,-1.); // flange
  volFlange -> SetLineColor(46); // rust

  volBypass
    = gGeoManager->MakeTube("FGBY",medMap.GetMedium(Geo::kCRGapBypass),
                             0,Geo::kFarBypassRad*scale,-1.); // bypass
  volBypass -> SetLineColor(kBlue);
  volBypass -> SetLineStyle(2); // dashed
  volBypass -> SetVisibility(kFALSE); // invisible by default

  volBypass -> AddNode(volFlange,1,gGeoIdentity);
  volFlange -> AddNode(volThroat,1,gGeoIdentity);
  volThroat -> AddNode(volCoil,1,new TGeoTranslation(0,-0.01*scale,0));

  return volBypass;

}

//_____________________________________________________________________________
TGeoVolume* GeoGeometry::BuildNearCoilAirGapVolume() const {
  //
  // Build near coil region to be used in airgaps between scint and steel
  // planes.  The near detector coil in the air gap
  // regions is built according to the same parameters as are used in
  // the build of the coil in the near detector full coverage scint pln:
  //   1)Insert a "bypass" tube (r=Geo::kNearFullBypassRad) region (not to 
  //     be confused with the bypass plastic insert).  This is constructed of 
  //     magnetized air.
  //   2)Insert into the "bypass" a "flange" (box of halfx/y = Geo::kFlangeRad)
  //     of  magnetized iron.
  //   3)Insert into the "flange" a "throat" (box of halfx/y = Geo::kThroatRad)
  //     of  magnetized air.
  //   4)Insert into the "throat" a "coil segment" (box of halfx/y
  //     = Geo::kCoilRad) of magnetized aluminum.  The coil segment is set
  //   5)Insert into the "coil segment" tubes (r=Geo::kNearCoolRad) of
  //     magnetized water in an array of 6 columns x 8 rows.
  //
  // The thickness in z of all pieces is set to be -1 when defining the shapes,
  // and built as a TGeoVolumeMulti.  When placed as a node in the scint
  // plane, the volumemulti shape will take on the thickness of its parent.
  // This allows for variable thickness of the scint plane in a non-perfect
  // geometry.
  //

  UpdateGlobalManager();

  const GeoMediumMap& medMap = GetMediumMap();
  
  // The nomenclature used here is "N" (near), "G" (gap), followed by
  // the two letter code to indicate the coil part type.
  std::string preface = "NG";

  // First check to see if volumemulti has been built. Need only build once.
  TGeoVolume* volBypass
              = gGeoManager->GetVolume(std::string(preface+"BY").c_str());
  if ( volBypass ) return volBypass;

  Double_t scale = Geo::GetScale(GetAppType());

  TGeoVolume* volWater
    = gGeoManager->MakeTube(std::string(preface+"WA").c_str(),
                            medMap.GetMedium(Geo::kCRGapWater),0,
                            Geo::kNearCoolRad*scale,-1.); // water
  volWater -> SetLineColor(kGreen);

  TGeoVolume* volCoil
    = gGeoManager->MakeBox(std::string(preface+"CO").c_str(),
                          medMap.GetMedium(Geo::kCRGapCoil),
                          Geo::kCoilRad*scale,Geo::kCoilRad*scale,-1.); // coil
  volCoil -> SetLineColor(kBlack);

  TGeoVolume* volThroat
    = gGeoManager->MakeBox(std::string(preface+"TR").c_str(),
                           medMap.GetMedium(Geo::kCRGapThroat),
                     Geo::kThroatRad*scale,Geo::kThroatRad*scale,-1.); //throat
  volThroat -> SetLineColor(kBlue);

  TGeoVolume* volFlange
    = gGeoManager->MakeBox(std::string(preface+"FL").c_str(),
                           medMap.GetMedium(Geo::kCRGapFlange),
                     Geo::kFlangeRad*scale,Geo::kFlangeRad*scale,-1.); //flange
  volFlange -> SetLineColor(kBlack);

  volBypass
      = gGeoManager->MakeTube(std::string(preface+"BY").c_str(),
                              medMap.GetMedium(Geo::kCRGapBypass),0,
                              Geo::kNearFullBypassRad*scale,-1.); // bypass
  volBypass -> SetLineColor(kBlue);
  volBypass -> SetLineStyle(2); // dashed

  volBypass -> AddNode(volFlange,1,gGeoIdentity);
  volFlange -> AddNode(volThroat,1,gGeoIdentity);

  volThroat -> AddNode(volCoil,1,
                       new TGeoTranslation(-0.01*scale,-0.01*scale,0));

  Int_t nx = 6;
  Int_t ny = 8;

  Float_t xedge = 0.02475*scale;// dist from x-edge to center of 1st(last) tube
  Float_t x0   = -Geo::kCoilRad*scale + xedge;  // first tube position in x
  // separation of tubes in x direction
  Float_t xsep = 2.*(Geo::kCoilRad*scale - xedge)/(Float_t)(nx-1);

  Float_t yedge = 0.02221*scale;// dist from y-edge to center of 1st(last) tube
  Float_t y0   = -Geo::kCoilRad*scale + yedge;  // first tube position in y
  // separation of tubes in y direction
  Float_t ysep = 2.*(Geo::kCoilRad*scale - yedge)/(Float_t)(ny-1);

  for ( int ix = 0; ix < nx; ix++ ) {
    Float_t xpos = x0 + (Float_t)ix*xsep;
    for ( int iy = 0; iy < ny; iy++ ) {
      Float_t ypos = y0 + (Float_t)iy*ysep;
      Int_t itube = ix*ny+iy+1;
      volCoil -> AddNode(volWater,itube,new TGeoTranslation(xpos,ypos,0));
    }
  }

  return volBypass;

}

---