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>
using namespace std;

#include "TGeoManager.h"
#include "TGeoNode.h"
#include "TGeoBBox.h"

#include "Conventions/Munits.h"
#include "GeoGeometry/Geo.h"
#include "GeoGeometry/GeoGeometry.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 "UgliGeometry/UgliDbiTables.h"
#include "UgliGeometry/UgliDbiStructHash.h"
#include "MessageService/MsgService.h"
#include "Plex/PlexVetoShieldHack.h"

ClassImp(GeoGeometry)

CVSID("$Id: GeoGeometry.cxx,v 1.22 2005/06/27 16:47:06 schubert Exp $");

//_____________________________________________________________________________
GeoGeometry::GeoGeometry() : fVldRange(),fAppType(Geo::kRecons),
                             fScale(1),fGeoManager(0) {
  // Default constructor

  UpdateGlobalManager();

  for ( int ism = 0; ism < 2; ism++ ) { 
     fSMZMin[ism] = +999.*fScale; 
     fSMZMax[ism] = -999.*fScale; 
  }
  
  MSG("Geo",Msg::kVerbose) << "GeoGeometry def ctor @ " << this << endl;

}

//_____________________________________________________________________________
GeoGeometry::GeoGeometry(const VldContext &vldc, Geo::EAppType apptype) : 
  fVldRange(),fAppType(apptype),fScale(1),fGeoManager(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.
  
  UpdateGlobalManager(); // set gGeoManager to null before building new 

  MSG("Geo",Msg::kVerbose) << "GeoGeometry normal ctor @ " << this << endl;
  
  //MsgStream* msgeo = MsgService::Instance() -> GetStream("Geo");
  //msgeo -> SetLogLevel(Msg::kDebug);
  
  for ( int ism = 0; ism < 2; ism++ ) { 
     fSMZMin[ism] = +999.*fScale; 
     fSMZMax[ism] = -999.*fScale; 
  }

  fScale = Geo::GetScale(apptype);

  TString name = vldc.AsString("s1ac-");
  fGeoManager = new TGeoManager(name.Data(),"a MINOS geometry");
  UpdateGlobalManager(); // set gGeoManager to fGeoManager
  
  // Build all materials and detector geometry
  this -> BuildAll(vldc);
  fGeoManager->CloseGeometry();

  TGeoVolume* volume = fGeoManager->GetTopVolume();
  UpdateNodeMatrices(volume);
  
  //this -> DumpVolume("MARS");
  
  fGeoManager->SetVisLevel(3);
  fGeoManager->SetVisOption(0);


}

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

  UpdateGlobalManager();

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

  fPlaneMap.clear();
  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 -> 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();
  
  TGeoVolume* linrVol = fGeoManager->GetVolume("LINR");
  GeoNode* hallNode = dynamic_cast<GeoNode*>(linrVol->GetNode(0));
  TGeoVolume* hallVol = hallNode -> GetVolume();
  TGeoBBox* hallBox = dynamic_cast<TGeoBBox*>(hallVol->GetShape());
  
  TVector3 localXYZ(-(hallBox->GetDX()),-(hallBox->GetDY()),
                    -(hallBox->GetDZ()));
  
  return (hallNode -> LocalToGlobal(localXYZ));

}

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

  UpdateGlobalManager();
  
  TGeoVolume* linrVol = fGeoManager->GetVolume("LINR");
  GeoNode* hallNode = dynamic_cast<GeoNode*>(linrVol->GetNode(0));
  TGeoVolume* hallVol = hallNode -> GetVolume();
  TGeoBBox* hallBox = dynamic_cast<TGeoBBox*>(hallVol->GetShape());
  
  TVector3 localXYZ(+(hallBox->GetDX()),+(hallBox->GetDY()),
                    +(hallBox->GetDZ()));
  return (hallNode -> LocalToGlobal(localXYZ));

}

//_____________________________________________________________________________
GeoSteelPlnNode* GeoGeometry::GetNearestSteelPlnNode(Double_t z) {
  // Public. Purpose:: Retrieve the nearest steel node for a given z
  // Notes: Method not const because GetSteelNodeMap builds map and stores as
  // data member.

  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;
  GeoSteelPlnNode* hiplnNode = 0;
  if ( steelplnItr != stnodemap.end() ) {
    hiplnNode = steelplnItr->second;
  }
  if ( steelplnItr != stnodemap.begin() ) {
    steelplnItr--;
    loplnNode = steelplnItr->second;
  }
  if ( hiplnNode == 0 ) return loplnNode;
  else if ( loplnNode == 0 ) return hiplnNode;
  
  // Bracketed by two planes, determine which is closer
  Double_t xyz_global[3] = {0,0,z};
  Double_t lo_xyz_local[3] = {0};       
  loplnNode->GlobalToLocal(xyz_global,lo_xyz_local);
  Double_t hi_xyz_local[3] = {0};       
  hiplnNode->GlobalToLocal(xyz_global,hi_xyz_local);

  if ( TMath::Abs(lo_xyz_local[2]) < TMath::Abs(hi_xyz_local[2]) ) {
    return loplnNode;
  }
  
  return hiplnNode;

}

//_____________________________________________________________________________
PlexPlaneId GeoGeometry::GetPlaneIdFromZ(Double_t z) const {
  // Purpose:: Retrieve detector plane corresponding to input z
  // This routine checks bounds of plane to determine if z is actually in
  // plane.
  // Use PlexPlaneId::IsValid() to determine validity of returned PlexPlaneId.
  // If returned plex is invalid => z not in plane.
 
  UpdateGlobalManager();
  
  MSG("Geo",Msg::kDebug) << "GeoGeometry::GetPlaneIdFromZ " << z << "."<< endl;
  assert(fGeoManager);

  PlexPlaneId invalidplex;

  // Move through the plane map from low to high z
  PlaneMapConstItr plnItr;
  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 ( TMath::Abs(xyz_local[2]) < dz ) {
      return plnItr->first;
    }
    else if ( xyz_local[2]  < -dz ) {
      // No possibility for intersection left assuming ordered planes from
      // lo to hi z
      return invalidplex;
    }
  }

  return invalidplex;

}

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

  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

  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() {
  // Purpose:: Return collection of ptrs for all scint pln nodes in detector

  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

  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() {
  // Public. Purpose:: Return collection of ptrs for all steel pln nodes 
  // in detector

  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() {
  // Public. Purpose:: Return collection of ptrs for all steel pln nodes 
  // in main body of detector, key'ed by z-position. 

  UpdateGlobalManager();

  if ( fSteelPlnNodeMap.empty() ) {
    PlaneMapConstItr plnItr;
    for ( plnItr = fPlaneMap.begin(); plnItr != fPlaneMap.end(); ++plnItr ){
      if ( (plnItr->first).IsSteel() ) {
        // For caldet, need to remove floor planes?
        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
 
  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

  UpdateGlobalManager();

  tmin = -999; tmax=-999;

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

//_____________________________________________________________________________
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

  UpdateGlobalManager();

  fGeoManager->ls(option);

}

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

  UpdateGlobalManager();

  std::string stoption(option);
  DumpVolume(stoption);
  
  cout << "GeoGeometry::Print\n" 
       << fVldRange.AsString() << endl;
  if ( fGeoManager ) 
    cout << "GeoManager " << fGeoManager->GetName() 
         << "," << fGeoManager->GetTitle() << endl;
  else 
    cout << "GeoManager not initialized." << endl;

}


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

  UpdateGlobalManager();

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

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

  if (vldc.GetDetector() == DetectorType::kUnknown) {
    MSG("Geo",Msg::kWarning)
      << "no possibility of building a geometry for " << endl
      << "   VldContext: " << vldc << endl;
    return;
  }

  // Build Detector Materials
  BuildMaterials(vldc);

  // Build Detector Geometry
  BuildGeometry(vldc);

}

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

  UpdateGlobalManager();

  MSG("Geo",Msg::kDebug) << "GeoGeometry BuildVldRange " << endl;
  
  // Start the VldRange covering all of time and then trim it down
  VldTimeStamp starttime = VldTimeStamp((time_t)0,0);
  VldTimeStamp endtime = VldTimeStamp(2038,1,18,0,0,0);
  fVldRange = VldRange(vldc.GetDetector(),vldc.GetSimFlag(),
                       starttime,endtime,"DBI");

  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());
  
}

//_____________________________________________________________________________
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::BuildMaterials(const VldContext& /* vldc */) {
   // Private. build this geometry's materials

   UpdateGlobalManager();

   MSG("Geo",Msg::kVerbose) << "GeoGeometry BuildMaterials " << endl;
   assert(fGeoManager);

//-----------List of Materials, Mixtures and Tracking Media--------------
   // Arguments: medName,medId,matId,
   //            isVol  (Sensitive volume flag),
   //            iField (User defined magnetic field, set = 2),
   //            fieldM (Max field value (in kGauss)),
   //            tMaxFd (Max angle due to field deflection(deg/step)),
   //            steMax (Max step allowed),
   //            deeMax (Max fractional energy loss in a step),
   //            epsil  (Tracking precision (cm)),
   //            stMin  (Minimum step sude to continuous processes (cm))

   TGeoMaterial *mat1 = new TGeoMaterial("ALUMINIUM",26.98,13,2.7);
   mat1->SetUniqueID(   1);
   new TGeoMedium("ALUMINIUM",1,1,0,0,0,20,0.1000000E+11,0.1829598,
                                       0.1000000E-03,0.2236820E-01);

   TGeoMaterial *mat2 = new TGeoMaterial("IRON",55.85,26,7.87);
   mat2->SetUniqueID(  2);
   new TGeoMedium("IRON",2,2,0,1,20,20,0.1000000E+11,0.25,
                                       0.1000000E-02,0.1952315E-01);
   new TGeoMedium("STEEL",3,2,0,0,20,20,0.1000000E+11,0.25,
                                       0.1000000E-02,0.1952315E-01);

   TGeoMixture *mat3 = new TGeoMixture("AIR",2,0.1205000E-02);
   mat3->SetUniqueID(  3);
   mat3->DefineElement(0,14.01,7,0.7);
   mat3->DefineElement(1,16.00,8,0.3);
   new TGeoMedium("AIR",4,3,0,0,0,20,0.1000000E+11,0.2488534,
                                          0.1000000E-03,0.7712014);

   TGeoMixture *mat4 = new TGeoMixture("CONCRETE",6,   2.50000    );
   mat4->SetUniqueID(  4);
   mat4->DefineElement(0,16,8,0.53);
   mat4->DefineElement(1,28.09,14,0.33);
   mat4->DefineElement(2,40.078,20,0.6300000E-01);
   mat4->DefineElement(3,22.99,11,0.1500000E-01);
   mat4->DefineElement(4,55.85,26,0.2000000E-01);
   mat4->DefineElement(5,26.98,13,0.4200000E-01);
   new TGeoMedium("CONCRETE",5,4,0,0,0,20,0.1000000E+11,0.1880915,
                                       0.1000000E-03,0.2108371E-01);

   TGeoMixture *mat5 = new TGeoMixture("ROCK",6,   2.50000    );
   mat5->SetUniqueID(  5);
   mat5->DefineElement(0,16,8,0.53);
   mat5->DefineElement(1,28.09,14,0.33);
   mat5->DefineElement(2,40.078,20,0.6300000E-01);
   mat5->DefineElement(3,22.99,11,0.1500000E-01);
   mat5->DefineElement(4,55.85,26,0.2000000E-01);
   mat5->DefineElement(5,26.98,13,0.4200000E-01);
   new TGeoMedium("ROCK",6,5,0,0,0,20,0.1000000E+11,0.1880915,
                                      0.1000000E-03,0.2108371E-01);

   TGeoMixture *mat6 = new TGeoMixture("PSTYRENE SCINT",2,   1.03200    );
   mat6->SetUniqueID(  6);
   mat6->DefineElement(0,1.00794,1,0.7742105E-01);
   mat6->DefineElement(1,12.011,6,0.9225789);
   new TGeoMedium("PSTYRENE SCINT",7,6,2,0,0,20,0.1000000E+11,0.3084481E-01,
                                                 0.1000000E-03,0.8242410E-02);

   TGeoMixture *mat7 = new TGeoMixture("COEX TIO2 PSTYRENE",4,   1.03200    );
   mat7->SetUniqueID( 7);
   mat7->DefineElement(0,1.00794,1,0.6967895E-01);
   mat7->DefineElement(1,12.011,6,0.8303211);
   mat7->DefineElement(2,44.88,22,0.5837760E-01);
   mat7->DefineElement(3,15.9994,8,0.4162240E-01);
   new TGeoMedium("COEX TIO2 PSTYRENE",8,7,0,0,0,20,0.1000000E+11,0.217354,
                                               0.1000000E-03,0.2307702E-01);

   TGeoMixture *mat8 = new TGeoMixture("FARCOIL",5,   2.34881    );
   mat8->SetUniqueID( 8);
   mat8->DefineElement(0,1.00794,1,0.1281415E-01);
   mat8->DefineElement(1,12.0107,6,0.6021570E-01);
   mat8->DefineElement(2,14.00674,7,0.1974151E-03);
   mat8->DefineElement(3,15.9994,8,0.2134690E-01);
   mat8->DefineElement(4,63.546,29,0.9054258);
   new TGeoMedium("FARCOIL",9,8,0,0,0,20,0.1000000E+11,0.1673649,
                                               0.1000000E-03,0.3520523E-01);


   return;

}

//_____________________________________________________________________________
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 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();

  // Pad x&y-dimensions because otherwise shield is cut-off 
  const Double_t xypad = 0.*fScale; // 0.2 meters
  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]) + xypad;
  Double_t yhsiz = 0.5*(hmaxsc[1] - hminsc[1]) + xypad;
  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
  TGeoMedium* medAIR    = fGeoManager->GetMedium("AIR");
  TGeoVolume* volHALL   = fGeoManager->MakeBox("HALL",medAIR,
                                                xhsiz,yhsiz,zhsiz);
  volHALL-> SetVisibility(kTRUE);
  volHALL-> SetLineColor(kGreen);

  // Concrete liner pad is box with halfwidth dimensions equal
  // to the hall halfwidths + the liner thickness (1 m)
  const Double_t linerpad = 1.0*fScale; // thickness concrete liner (1 m)
  TGeoMedium* medCONCRETE = fGeoManager->GetMedium("CONCRETE");
  TGeoVolume* volLINR     = fGeoManager->MakeBox("LINR",medCONCRETE,
                                         xhsiz+linerpad,yhsiz+linerpad,
                                         zhsiz+linerpad);
  volLINR-> SetVisibility(kTRUE);
  volLINR-> SetLineColor(kYellow);

  // Rock shell
  // rock padding is 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]-xypad),
                                TMath::Abs(hmaxsc[0]+xypad));
  Double_t absymax = TMath::Max(TMath::Abs(hminsc[1]-xypad),
                                TMath::Abs(hmaxsc[1]+xypad));
  Double_t abszmax = TMath::Max(TMath::Abs(hminsc[2]),
                                TMath::Abs(hmaxsc[2]));
  TGeoMedium* medROCK = fGeoManager->GetMedium("ROCK");
  TGeoVolume* volMARS = fGeoManager->MakeBox("MARS",medROCK,
                                              absxmax+rockpad,absymax+rockpad,
                                              abszmax+rockpad);
  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;
  new GeoNode(this,volLINR,new TGeoTranslation("LINR",x0hall,y0hall,z0hall),
              volMARS,globalpath,nodename);
  
  nodename = "HALL_1";
  globalpath += "/" + nodename;
  GeoNode* hallNode = new GeoNode(this,volHALL,gGeoIdentity,
                                  volLINR,globalpath,nodename);

  BuildDetector(vldc,hallNode);
  
}

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

  UpdateGlobalManager();

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

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

  TGeoVolume* hallVol  = hallNode -> GetVolume();
  std::string hallPath = hallNode -> GetGlobalPath();
  
  // 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(0);
  const Double_t* h0 = linrNode -> GetMatrix() -> GetTranslation();
  
  TGeoMedium* steelMed = fGeoManager -> GetMedium("IRON"); // has b-field
  TGeoMedium* airMed   = fGeoManager -> GetMedium("AIR");

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

  Int_t plnMin[2] = {-1};
  Int_t plnMax[2] = {-1};
  Int_t nSM = GetSMPlaneLimits(vldc,plnMin,plnMax);

  
  // detx0, dety0 are used to position the coil
  Float_t detx0 = 0;
  Float_t dety0 = 0;

  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() ) 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;
    
    if ( irow % 20 == 0 ) {
      MSG("Geo",Msg::kInfo) << "Building steel row "<< stName.c_str()<< endl;
    } 
           
    Float_t pairhalfx = 0;
    Float_t pairhalfy = 0;
    Float_t pairhalfz = 0.5*(stRow->GetTotalZ())*fScale;

    if ( !steelId.IsVetoShield() ) {
      stplnVolume = new GeoSteelPlnVolume(this,steelId,
                                        stRow->GetThickness()*fScale,steelMed);
      stplnVolume -> SetLineColor(kRed);
      stplnVolume -> SetVisibility(kTRUE);

      pairhalfx = dynamic_cast<TGeoBBox*>(stplnVolume->GetShape())->GetDX();
      // This fudge is because bbox of neardetector planes only includes
      // plane width and not near detector offset which seems worrisome
      // given that this offset (asymmetry about x) is included in the
      // shape definition.
      if ( vldc.GetDetector() == Detector::kNear ) pairhalfx += 
                                        fScale*TMath::Abs(Geo::kNearXOffset);
      pairhalfy = dynamic_cast<TGeoBBox*>(stplnVolume->GetShape())->GetDY();
    }
    
    // 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() ) {
      MSG("Geo",Msg::kVerbose) << "scint pln name " << scName.c_str() 
                               << endl;
      scRow = ugliDbiTables.GetDbiScintPlnByIndex(scintId.GetPlane());
    
      scplnVolume = new GeoScintPlnVolume(this,scintId,
                                          scRow->GetThickness()*fScale,
                                          airMed);
      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.
      pairhalfx = TMath::Max(pairhalfx,(Float_t)dynamic_cast<TGeoBBox*>
                             (scplnVolume->GetShape())->GetDX());
      pairhalfy = TMath::Max(pairhalfy,(Float_t)dynamic_cast<TGeoBBox*>
                             (scplnVolume->GetShape())->GetDY());
    }

    // Build pair bounding box and position according to steel plane coords
    std::string pairName = GetGeoCompatibleName(steelId.AsString("b"));
    TGeoVolume* pairVol = fGeoManager->MakeBox(pairName.c_str(),airMed,
                                               pairhalfx,pairhalfy,pairhalfz);
    pairVol -> SetVisibility(kFALSE);
    pairVol -> VisibleDaughters(kTRUE);

    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 and not hall, move to hall 
    // local
    Float_t x0pair = (stRow->GetX0())*fScale - h0[0];
    Float_t y0pair = (stRow->GetY0())*fScale - h0[1];
    Float_t z0pair = (stRow->GetZBack()-0.5*(stRow->GetTotalZ()))*fScale-h0[2];
    TGeoCombiTrans* pairMatrix = new TGeoCombiTrans(pairName.c_str(),
                                             x0pair,y0pair,z0pair,stRotMatrix);
    MSG("Geo",Msg::kDebug) << "Placing plane pair at " << x0pair << ","
                             << y0pair << "," << z0pair << endl;

    std::string pairPath = hallPath + "/" + pairName+"_1";
    new GeoNode(this,pairVol,pairMatrix,hallVol,pairPath,pairName+"_1");

    if ( stplnVolume ) {
      // Place the steel plane within the pairVol
      std::string stplnNodeName = stName + "_1";
      std::string steelPlnPath = pairPath + "/" + stplnNodeName;
      new GeoSteelPlnNode(this,stplnVolume,
                          new TGeoTranslation(stName.c_str(),0.,0.,
          0.5*(stRow->GetTotalZ()-stRow->GetThickness())*fScale),pairVol,
                          steelPlnPath,stplnNodeName,steelId);
      int iSM = 0;
      if ( nSM > 1 && steelId.GetPlane() > plnMax[0] ) iSM = 1;
      detx0 = x0pair; // hall coordinates
      dety0 = y0pair; // hall coordinates
      // fSMZMin/Max in MARS coordinates
      fSMZMin[iSM] = TMath::Min(fSMZMin[iSM],(stRow->GetZBack()
                              -stRow->GetTotalZ())*fScale);  
      fSMZMax[iSM] = TMath::Max(fSMZMax[iSM],(stRow->GetZBack())*fScale);
      
    }
    
    if ( scplnVolume ) {
      // Place the scint plane within the pairVol
      Float_t x0pln = (scRow->GetX0RelSteel())*fScale;
      Float_t y0pln = (scRow->GetY0RelSteel())*fScale;
      Float_t z0pln = 0.5*(-stRow->GetTotalZ()+scRow->GetThickness())*fScale;
      
      MSG("Geo",Msg::kVerbose) << "Placing scint pln at " << x0pln << ", "
                               << y0pln << "," << z0pln << endl;
        
      // Rotation matrix relative to steel
      // (thetax,phix) = (90, 0+zrotrelsteel)
      // (thetay,phiy) = (90,90+zrotrelsteel)
      // (thetaz,phiz) = ( 0, 0)
      Double_t zrotrelsteel = scRow->GetZRotRelSteelDeg();
      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);
    } 
  } // end of loop over all plane db rows

  // Build coil after all planes are positioned so that detector
  // limits are known
  DetectorType::Detector_t dettype = (DetectorType::Detector_t)
                                      vldc.GetDetector();
  switch ( dettype ) {
  case DetectorType::kFar:
    BuildFarCoil(hallNode,detx0,dety0);
    break;
  case DetectorType::kNear:
    BuildNearCoil(hallNode,detx0,dety0);
    break;
  default:
    break;
  }

  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.

  DetectorType::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"
    << DetectorType::AsString(dettype) << "." << endl;
  } // end of switch

  return nSM;
  
}

//_____________________________________________________________________________
void GeoGeometry::BuildNearCoil(const TGeoNode* hallNode,
                                Float_t x0, Float_t y0) {
  // 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);

  TGeoMedium* medSTEEL    = fGeoManager->GetMedium("STEEL");
  TGeoMedium* medAIR      = fGeoManager->GetMedium("AIR");
  TGeoMedium* medAL       = fGeoManager->GetMedium("ALUMINIUM");

  TGeoVolume* hallVol = hallNode -> GetVolume();

  // The dimensions of the coil should be extracted from the db
  Float_t shalfx  = Geo::kNearSteelHoleHalfWidth*fScale; // steel, 0.18 m 
  Float_t shalfy  = shalfx; // steel
  Float_t ahalfx  = 0.13*fScale; // air, halfwidth 0.13 m
  Float_t ahalfy  = ahalfx; // air
  Float_t alhalfx = 0.1118*fScale; // aluminum
  Float_t alhalfy = 0.1143*fScale; // 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(0);
  const Double_t* h0 = linrNode -> GetMatrix() -> GetTranslation();
  zmin -= h0[2];
  zmax -= h0[2];
  
  Float_t zcoillen = zmax - zmin;
  TGeoVolume* BUS0 = fGeoManager->MakeBox("BUS0",medSTEEL,
                                           shalfx,shalfy,0.5*zcoillen);
  BUS0 -> SetLineColor(14);  // gray  
  BUS0 -> SetVisibility(kTRUE);

  TGeoVolume* BUA0 = fGeoManager->MakeBox("BUA0",medAIR,
                                           ahalfx,ahalfy,0.5*zcoillen);
  BUA0 -> SetLineColor(14);  
  BUA0 -> SetVisibility(kTRUE);
  BUS0 -> AddNode(BUA0,1,gGeoIdentity);

  TGeoVolume* BUC0 = fGeoManager->MakeBox("BUC0",medAL,
                                           alhalfx,alhalfy,0.5*zcoillen);
  BUC0 -> SetLineColor(14);  
  BUC0 -> SetVisibility(kTRUE);
  BUA0 -> AddNode(BUC0,1,gGeoIdentity);

  // Place inner coil in hall
  // Rotate first, then translate
  TGeoRotation *rot45 = new TGeoRotation("BUS0",90,45,90,135,0,0);
  rot45 -> RegisterYourself();
  Float_t xcoil = x0;
  Float_t ycoil = y0;
  Float_t zcoil = (zmin+zmax)/2.;
  hallVol -> AddNode(BUS0,1,new TGeoCombiTrans(xcoil,ycoil,zcoil,rot45));

  // Place return coil in hall
  Float_t steelhalfy = 1.9*fScale; // 1.9 m halfwidth in y direction
  xcoil = x0 - steelhalfy;
  ycoil = y0 - steelhalfy;
  hallVol -> AddNode(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",medSTEEL,
                                            shalfx,0.5*coillenxy,shalfy);
  BUSXY -> SetLineColor(14);  // gray  
  BUSXY -> SetVisibility(kTRUE);
  TGeoVolume* BUAXY = fGeoManager->MakeBox("BUAXY",medAIR,
                                            ahalfx,0.5*coillenxy,ahalfy);
  BUAXY -> SetLineColor(14);  
  BUAXY -> SetVisibility(kTRUE);
  BUSXY -> AddNode(BUAXY,1,gGeoIdentity);

  TGeoVolume* BUCXY = fGeoManager->MakeBox("BUCXY",medAL,
                                            alhalfx,0.5*coillenxy,alhalfy);
  BUCXY -> SetLineColor(14);  
  BUCXY -> SetVisibility(kTRUE);
  BUAXY -> AddNode(BUCXY,1,gGeoIdentity);

  // Will be rotated and then translated
  TGeoRotation *rotneg45 = new TGeoRotation("BUSXY",90,315,90,45,0,0);
  rotneg45 -> RegisterYourself();
  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(const TGeoNode* hallNode, 
                               Float_t x0, Float_t y0) {
  // Private. Build the far detector coil. Fix me.  
  // x0,y0 is coil hole center in HALL coordinates (assumed constant
  // over all planes).

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

  // Guess vertical drop of coil along end planes - from y-center of coil hole
  // to y-center of return coil on floor.
  Float_t deltay = 5.*fScale; // 5 m drop, guess

  // These dimensions/etc. of the coil should be in the db, but for now
  // use hardwired values
  TGeoMedium* medFARCOIL = fGeoManager->GetMedium("FARCOIL");
  TGeoMedium* medSTEEL   = fGeoManager->GetMedium("STEEL");

  TGeoVolume* hallVol = hallNode -> GetVolume();

  // Steel shell surrounding copper coil 
  Float_t sRad = Geo::kFarSteelHoleRad*fScale;  // 15 cm
  Float_t cRad = 0.125*fScale; // copper coil radius, 12.5 cm, tdr

  // SM0 coil along z-direction of detector, ends at supermodule ends
  Float_t zmin0,zmax0;
  GetZExtent(zmin0,zmax0,0);
  // Convert to HALL coordinates
  TGeoNode* linrNode = fGeoManager->GetVolume("MARS")->GetNode(0);
  const Double_t* h0 = linrNode -> GetMatrix() -> GetTranslation();
  zmin0 -= h0[2];
  zmax0 -= h0[2];

  Float_t zcoillen0 = zmax0-zmin0;
  TGeoVolume* BUS0=fGeoManager->MakeTube("BUS0",medSTEEL,0,sRad,0.5*zcoillen0);
  BUS0 -> SetLineColor(14);  // gray  
  BUS0 -> SetVisibility(kTRUE);
  TGeoVolume* BUC0=fGeoManager->MakeTube("BUC0",medFARCOIL,0,cRad,
                                         0.5*zcoillen0);
  BUC0 -> SetLineColor(14);  
  BUC0 -> SetVisibility(kTRUE);
  BUS0 -> AddNode(BUC0,1,gGeoIdentity);

  // Place SM0 inner coil in hall
  Float_t z0 = 0.5*(zmin0 + zmax0);
  hallVol -> AddNode(BUS0,1,new TGeoTranslation("BUS0",x0,y0,z0));
  // Place SM0 return coil in hall
  hallVol -> AddNode(BUC0,1,new TGeoTranslation("BUC0",x0,y0-deltay,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->MakeTube("BUS1",medSTEEL,0,sRad,0.5*zcoillen1);
  BUS1 -> SetLineColor(14);  // gray
  BUS1 -> SetVisibility(kTRUE);
  TGeoVolume* BUC1=fGeoManager->MakeTube("BUC1",medFARCOIL,
                                          0,cRad,0.5*zcoillen1);
  BUC1 -> SetLineColor(14);
  BUC1 -> SetVisibility(kTRUE);
  BUS1 -> AddNode(BUC1,1,gGeoIdentity);

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

  // Build coil volume on xy faces
  Float_t ycoillen = deltay + 2.*cRad;
  TGeoVolume* BUXY=fGeoManager->MakeTube("BUXY",medFARCOIL,
                                          0,cRad,0.5*ycoillen);
  BUXY -> SetLineColor(14); // gray  
  BUXY -> SetVisibility(kTRUE);

  // Will be rotated and then translated
  TGeoRotation *rot90 = new TGeoRotation("BUXY",90,0,180,0,90,90);
  rot90 -> RegisterYourself();

  Float_t xcoil = x0;
  Float_t ycoil = (y0 + cRad) - ycoillen/2.;
  hallVol -> AddNode(BUXY,1,
                  new TGeoCombiTrans("BUXY_1",xcoil,ycoil,zmin0-cRad,rot90));
  hallVol -> AddNode(BUXY,2,
                  new TGeoCombiTrans("BUXY_2",xcoil,ycoil,zmax0+cRad,rot90));
  hallVol -> AddNode(BUXY,3,
                  new TGeoCombiTrans("BUXY_3",xcoil,ycoil,zmin1-cRad,rot90));
  hallVol -> AddNode(BUXY,4,
                  new TGeoCombiTrans("BUXY_4",xcoil,ycoil,zmax1+cRad,rot90));

  return;

}

//_____________________________________________________________________________
void GeoGeometry::BuildModules(const UgliDbiTables& ugliDbiTables,
                               GeoScintPlnNode* plnNode) {
  // Private. Purpose: Build the modules as part of the scintillator plane
  
  UpdateGlobalManager();
  
  assert(fGeoManager);
  
  TGeoMedium* medAL  = fGeoManager -> GetMedium("ALUMINIUM");
  TGeoMedium* medAIR = fGeoManager -> GetMedium("AIR");

  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,medAL,
                                                      medAIR);
 
    // 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 = stpName+"_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);

  TGeoMedium* scintMed = fGeoManager -> GetMedium("PSTYRENE SCINT");
  TGeoMedium* coexMed  = fGeoManager -> GetMedium("COEX TIO2 PSTYRENE");

  DbiResultPtr<UgliDbiStripStruct> stripTbl(vldc);
  
  for ( unsigned int irow = 0; irow < stripTbl.GetNumRows(); irow++ ) {
    const UgliDbiStripStruct* sRow = stripTbl.GetRow(irow);
    string stpName = GetGeoCompatibleName(sRow -> GetShapeName());

    // 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,
                                        scintMed,coexMed);
  } 

}

//_____________________________________________________________________________
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);
  }

}

---