GeoScintMdlVolume.cxx

Go to the documentation of this file.

// $Id: GeoScintMdlVolume.cxx,v 1.18 2008/01/09 06:38:11 schubert Exp $
//
// GeoScintMdlVolume
//
// GeoScintMdlVolume is a scint module volume. The volume definition
// consists of an aluminum shell with an air interior.
//
// Author:  S. Kasahara 07/04

#include <iostream>
using namespace std;
#include <TGeoManager.h>
#include <TGeoMatrix.h>
#include <TGeoBBox.h>
#include <TGeoArb8.h> // has TGeoTrap definition
#include <TGeoXtru.h>
#if ROOT_VERSION_CODE >= ROOT_VERSION(5,6,0)
#include <TGeoScaledShape.h>
#endif
#include <TROOT.h>
#include <TMath.h>

#include "GeoGeometry/GeoScintMdlVolume.h"
#include "GeoGeometry/GeoGeometry.h"
#include "GeoGeometry/GeoMediumMap.h"
#include "GeoGeometry/GeoStripVolume.h"
#include "GeoGeometry/GeoStripNode.h"
#include "GeoGeometry/GeoMedium.h"
#include "MessageService/MsgService.h"
#include "UgliGeometry/UgliDbiTables.h"
#include "UgliGeometry/UgliDbiStructHash.h"

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

//_____________________________________________________________________________
GeoScintMdlVolume::GeoScintMdlVolume() : GeoVolume(),fScintMdlId(),fThick(0),
                                         fSkinThick(0),
                                         fBaseShp(0),fBaseShpAir(0) {

  if ( gROOT -> GetVersionInt() < 50600 ) {
    MSG("Geo",Msg::kFatal) << "GeoGeometry requires ROOT version >= 5.06/00."
                           << "\n Please switch to a newer version of root." 
                           << endl;
    abort();
  }  
    
  // Default constructor, for i/o
  for ( int ip = 0; ip < 4; ip++ ) {
    for ( int ic = 0; ic < 2; ic++ ) fXY[ip][ic] = 0;
  }

}

//_____________________________________________________________________________
GeoScintMdlVolume::GeoScintMdlVolume(GeoGeometry* geo,
                                     const PlexScintMdlId& mdlId,
                                     const UgliDbiTables& ugliTables) :
    GeoVolume(geo,
        (GeoGeometry::GetGeoCompatibleName(mdlId.AsString())).c_str(),0),
    fScintMdlId(mdlId),fThick(0),fSkinThick(0),fBaseShp(0),fBaseShpAir(0) {
  // Normal constructor - Protected.

  if ( gROOT -> GetVersionInt() < 50600 ) {
    MSG("Geo",Msg::kFatal) << "GeoGeometry requires ROOT version >= 5.06/00."
                           << "\n Please switch to a newer version of root." 
                           << endl;
    abort();
  }  

  UpdateGlobalManager();

  const GeoMediumMap& medMap = fGeoGeometry -> GetMediumMap();
  if ( fScintMdlId.IsVetoShield() ) 
                 SetMedium(medMap.GetMedium(Geo::kVSMdlAlSkin));
  else SetMedium(medMap.GetMedium(Geo::kMdlAlSkin));

  MSG("Geo",Msg::kDebug) << "GeoScintMdlVolume ctor, mdl name " 
                         << GetName() << endl;
 
  fShpName = GeoGeometry::GetGeoCompatibleName(mdlId.AsString());
   
  Float_t scale = GetScale();
  fSkinThick = scale*Geo::kModuleSkinThick;  // thickness of module skin
  fThick = scale*Geo::kStripThick + 2.*fSkinThick; // thickness of module
  
  // "s" requests shp name
  std::string baseshpname 
              = GeoGeometry::GetGeoCompatibleName(mdlId.AsString("s"));
  std::string baseshpairname = baseshpname + "_air";
  fBaseShp = dynamic_cast<TGeoShape*>(gGeoManager->GetListOfShapes()
                            -> FindObject(baseshpname.c_str()));
  fBaseShpAir = dynamic_cast<TGeoShape*>(gGeoManager->GetListOfShapes()
                            -> FindObject(baseshpairname.c_str()));

  if ( !fBaseShp ) {
    // create base shape if first call
    fBaseShp = BuildModuleShape(ugliTables,baseshpname);
    fBaseShpAir = BuildModuleShape(ugliTables,baseshpairname,true);
  }
  
  TGeoShape* shape = BuildModuleShape(ugliTables,fShpName);

  SetShape(shape);

  // Because shape was unknown at time of TGeoVolume construction, volume
  // was added to run-time list of volumes. Remove and add to normal list
  gGeoManager -> GetListOfGVolumes() -> Remove(this);
  Int_t index = gGeoManager->GetListOfVolumes() -> GetEntriesFast();
  gGeoManager -> GetListOfVolumes() -> AddAtAndExpand(this,index);
  
  this -> SetLineColor(17); // Grey
  //this -> SetLineWidth(GetLineWidth()*2);
  this -> SetVisibility(kTRUE);

  this -> AddAirNode(ugliTables);
  
}

//_____________________________________________________________________________
void GeoScintMdlVolume::AddAirNode(const UgliDbiTables& ugliTables) {  
  // Build shape of module interior and associated air volume. Place
  // air volume as node of aluminum module.

  UpdateGlobalManager();

  Geo::EDetComponent detcomp = Geo::kMdlAir;
  if ( fScintMdlId.IsVetoShield() ) detcomp = Geo::kVSMdlAir;

  std::string airShpName = fShpName+"_air";
  TGeoShape* airShape = BuildModuleShape(ugliTables,airShpName,true);
  
  std::string airVolName = std::string(GetName()) + "_air"; 
  TGeoVolume* airVol = new TGeoVolume(airVolName.c_str(),airShape,
                       fGeoGeometry->GetMediumMap().GetMedium(detcomp));
  airVol -> SetLineColor(17); // 17 for Grey
  airVol -> SetVisibility(kTRUE);
  
  this -> AddNode(airVol,1,gGeoIdentity);

  return;
                                    
}

//_____________________________________________________________________________
void GeoScintMdlVolume::Print(Option_t* /* option */) const {
   // print to cout

  UpdateGlobalManager();

  cout << "GeoScintMdlVolume::Print\n"
       << GetName() << endl;
  
}

//_____________________________________________________________________________
GeoStripNode* GeoScintMdlVolume::GetStripNode(const PlexStripEndId& seid) 
                                                                       const {
   // Retrieve strip node corresponding to this strip end id

  UpdateGlobalManager();

  // Strip sub-part and/or end info off seid for comparison
  PlexStripEndId seid_clean(seid);
  seid_clean.SetEnd(StripEnd::kUnknown);
  seid_clean.SetSubPart(StripEnd::kWhole);

  TGeoVolume* airVol = GetAirNode()->GetVolume();
  Int_t nstrip = airVol->GetNdaughters();
  for ( int istp = 0; istp < nstrip; istp++ ) {
    GeoStripNode* stpnode = dynamic_cast<GeoStripNode*>(airVol->GetNode(istp));
    if ( seid_clean == stpnode->GetSEId() ) {
      return stpnode;
    }
  }

  return 0;
  
}

//_____________________________________________________________________________
std::vector<GeoStripNode*> GeoScintMdlVolume::GetStripNodePtrVector() const {
   // Retrieve vector of strip nodes in this module

  UpdateGlobalManager();

  std::vector<GeoStripNode*> stpnodes;
  TGeoVolume* airVol = GetAirNode()->GetVolume();
  Int_t nstrip = airVol->GetNdaughters();
  for ( int istp = 0; istp < nstrip; istp++ ) {
    GeoStripNode* stpnode = dynamic_cast<GeoStripNode*>(airVol->GetNode(istp));
    stpnodes.push_back(stpnode);
  }

  return stpnodes;
  
}


//_____________________________________________________________________________
TGeoShape* GeoScintMdlVolume::BuildModuleShape(const UgliDbiTables& ugliTables,
                                               std::string shpName, 
                                               bool isInner) {
  // Private method used to build shape for this plane. 

  UpdateGlobalManager();
  
  MSG("Geo",Msg::kDebug) << "BuildModuleShape for shape " << shpName.c_str() 
                         << endl;
  
  
  TGeoShape* mdlBaseShp = fBaseShp;
  if ( isInner ) {
    mdlBaseShp = fBaseShpAir;
  }

  if ( mdlBaseShp ) {
    BuildHalfY(ugliTables); // simpler than Build4Corners
    Double_t halfybasic = dynamic_cast<TGeoBBox*>(mdlBaseShp) -> GetDY();
    Double_t halfy = 0.5*(fXY[1][1] - fXY[0][1]);
    if ( isInner ) halfy -= fSkinThick;
    Double_t scaley = halfy/halfybasic;
    if ( TMath::Abs(scaley - 1.) > 1.e-6 ) {
#if ROOT_VERSION_CODE >= ROOT_VERSION(5,6,0)
      return (new TGeoScaledShape(shpName.c_str(),mdlBaseShp,
                                  new TGeoScale(1.,scaley,1.)));
#endif
    }
    else {
      return mdlBaseShp;
    }
  }

  TGeoShape* shape = 0;
  
  Detector::Detector_t dettype = fScintMdlId.GetDetector();

  switch ( dettype ) {

  case Detector::kFar:

    shape = BuildFar(ugliTables,shpName,isInner);
    break;
    
  case Detector::kNear:

    shape = BuildNear(ugliTables,shpName,isInner);
    break;

  case Detector::kCalDet:

    shape = BuildCalDet(ugliTables,shpName,isInner);
    break;

  default:

    MSG("Geo",Msg::kError)<< "Scint mdl shape construction for detector type\n"
    << Detector::AsString(dettype)
    << " not yet supported." << endl;

  } // end of switch

  return shape;

} 

//_____________________________________________________________________________
TGeoShape* GeoScintMdlVolume::BuildCalDet(const UgliDbiTables& ugliTables,
                                          std::string shpName, bool isInner) {
  // Private method used to build shape for calibration detector modules. 

  UpdateGlobalManager();

  UChar_t mdlType = fScintMdlId.GetModuleType();
  std::string regionName = ( isInner ) ? "interior" : "shell";
  
  MSG("Geo",Msg::kVerbose)
    << "Build cal detector module shape for module type " 
    << mdlType << " " << regionName.c_str() << "." << endl;

  TGeoShape* shape = 0;

  switch ( mdlType ) {
  default:
    // All caldet module types are simple boxes
    shape = BuildBBox(ugliTables,shpName,isInner);
    break;
  } // end of switch

  return shape;

} 

//_____________________________________________________________________________
TGeoShape* GeoScintMdlVolume::BuildFar(const UgliDbiTables& ugliTables,
                                       std::string shpName, bool isInner) {
  // Private method used to build shape for far detector modules. 

  UpdateGlobalManager();

  UChar_t mdlType = fScintMdlId.GetModuleType();
  std::string regionName = ( isInner ) ? "interior" : "shell";
  
  MSG("Geo",Msg::kVerbose)
    << "Build far detector module shape for module type " 
    << mdlType << " " << regionName.c_str() << "." << endl;

  TGeoShape* shape = 0;

  switch ( mdlType ) {
    case 'A':
    case 'a':
    case 'B':
    case 'b':
      // Simple trapezoid
      shape = BuildSimpleTrap(ugliTables,shpName,isInner);
      break;
    case 'C':
    case 'c':
    case 'E':
    case 'e':
      // Box 
      shape = BuildBBox(ugliTables,shpName,isInner);
      break;
    case 'd':
    case 'f':
      // Box with coil hole cut out on high tpos side
      shape = BuildFardfXtru(ugliTables,shpName,isInner);
      break;
    case 'D':
    case 'F':
      // Box with coil hole cut out on lo tpos side
      shape = BuildFarDFXtru(ugliTables,shpName,isInner);
      break;
    default:
      MSG("Geo",Msg::kError)
      << "Scint module shape construction for far detector module type\n"
      << mdlType << " not yet supported." << endl;
  } // end of switch

  return shape;

} 

//_____________________________________________________________________________
TGeoShape* GeoScintMdlVolume::BuildNear(const UgliDbiTables& ugliTables,
                                        std::string shpName, bool isInner) {
  // Private method used to build shape for near detector modules. 

  UpdateGlobalManager();

  UChar_t mdlType = fScintMdlId.GetModuleType();
  std::string regionName = ( isInner ) ? "interior" : "shell";
  
  MSG("Geo",Msg::kVerbose)
    << "Build near detector module shape for module type " 
    << mdlType << " " << regionName.c_str() << "." << endl;

  TGeoShape* shape = 0;

  switch ( mdlType ) {
    case 'g':
    case 'G':
    case 'h':
    case 'H':
    case 'j':
    case 'J':
    case 'k':
    case 'K':
    case 'n':
    case 'N':
       // Trapezoids or box
      shape = BuildNearXtru(ugliTables,shpName,isInner);
      break;
    case 'i':
     // Box with coil hole in lower tpos, lower lpos corner
      shape = BuildNeariXtru(ugliTables,shpName,isInner);
      break;
    case 'I':
      // Box with coil hole in upper tpos, lower lpos corner
      shape = BuildNearIXtru(ugliTables,shpName,isInner);
      break;
    case 'l':
    case 'M':
      // Box or trapezoid with coil hole on lo tpos side, centered on lpos=0
      shape = BuildNearlMXtru(ugliTables,shpName,isInner);
      break;
    case 'L':
    case 'm':
      // Box or trapezoid with coil hole on hi tpos side, centered on lpos=0
     shape = BuildNearLmXtru(ugliTables,shpName,isInner);
     break;
    default:
      MSG("Geo",Msg::kError)
      << "Scint module shape construction for near detector module type\n"
      << mdlType << " not yet supported." << endl;
  } // end of switch

  return shape;

} 

//_____________________________________________________________________________
TGeoShape* GeoScintMdlVolume::BuildSimpleTrap(const UgliDbiTables& ugliTables,
                                              std::string shpName, 
                                              bool isInner) {
 // Simple Trapezoid shape

  UpdateGlobalManager();

  Float_t m01 =0; // slope on lo lpos side
  Float_t m23 = 0; // slope on hi lpos side
  Build4Corners(ugliTables,m01,m23); // stores 4 corners of trap in fXY

  Double_t halfxloy = 0.5*(fXY[3][0] - fXY[0][0]);
  Double_t halfxhiy = 0.5*(fXY[2][0] - fXY[1][0]);
  Double_t halfy = 0.5*(fXY[1][1] - fXY[0][1]);
  Double_t halfz = 0.5*fThick;

  if ( isInner )  {
    // Interior shape
    Double_t b01 = fXY[1][1] - m01*fXY[1][0];
    Double_t arctanm = atan(TMath::Abs(m01));
    Double_t delta = fSkinThick/sin(arctanm);
    Double_t y0p = fXY[0][1] + fSkinThick;
    Double_t x0p = (y0p - b01)/m01 + delta;
    Double_t y1p = fXY[1][1] - fSkinThick;
    Double_t x1p = (y1p - b01)/m01 + delta;
    
    halfxloy = -x0p;
    halfxhiy = -x1p;
    halfy -= fSkinThick;
    halfz -= fSkinThick;
  }

  TGeoTrap* mdlShape = new TGeoTrap(shpName.c_str(), halfz, 0., 0.,
                                    halfy,halfxloy,halfxhiy,0.,halfy,
                          halfxloy,halfxhiy,0.);
  
  return mdlShape;

} 

//_____________________________________________________________________________
TGeoShape* GeoScintMdlVolume::BuildBBox(const UgliDbiTables& ugliTables,
                                        std::string shpName, bool isInner) {
 // Simple Box shape

  UpdateGlobalManager();

  Float_t m01 =0; // slope on lo lpos side
  Float_t m23 = 0; // slope on hi lpos side
  Build4Corners(ugliTables,m01,m23); // stores 4 corners of trap in fXY

  Double_t halfx = 0.5*(fXY[3][0] - fXY[0][0]);
  Double_t halfy = 0.5*(fXY[1][1] - fXY[0][1]);
  Double_t halfz = 0.5*fThick;

  if ( isInner )  {
    // Interior shape
    halfx -= fSkinThick;
    halfy -= fSkinThick;
    halfz -= fSkinThick;
  }

  TGeoBBox* mdlShape = new TGeoBBox(shpName.c_str(), halfx, halfy, halfz);
  
  return mdlShape;

} 

//_____________________________________________________________________________
TGeoShape* GeoScintMdlVolume::BuildFardfXtru(const UgliDbiTables& ugliTables,
                                             std::string shpName,
                                             bool isInner) {
 // Box with coil hole cut out from top

  UpdateGlobalManager();

  Float_t m01 =0; // slope on lo lpos side
  Float_t m23 = 0; // slope on hi lpos side
  Build4Corners(ugliTables,m01,m23); // stores 4 corners of trap in fXY

  const int np = 15; // greater than 5
  Double_t xp[np] = {0}, yp[np] = {0};

  Double_t halfthick = 0.5*fThick;
  Double_t skinthick = 0;
  Double_t radhole = Geo::kFarModuleHoleRad*GetScale();
  Double_t delta = 0;
  Double_t dang = TMath::Pi()/(np-5);
  radhole /= cos(0.5*dang);
  
  if ( isInner ) {
    halfthick -= fSkinThick;
    skinthick = fSkinThick;
    radhole += fSkinThick;
    delta = TMath::ASin(skinthick/radhole); // small angle approx
    dang = (TMath::Pi()-2.*delta)/(np-5);
  }
    
  xp[0] = fXY[0][0] + skinthick; yp[0] = fXY[0][1] + skinthick;
  xp[1] = fXY[1][0] + skinthick; yp[1] = fXY[1][1] - skinthick;
  xp[np-2] = fXY[2][0] - skinthick; yp[np-2] = fXY[2][1] - skinthick;
  xp[np-1] = fXY[3][0] - skinthick; yp[np-1] = fXY[3][1] + skinthick;
  Double_t holevtx = fXY[2][1];
  for ( int ip = 2; ip < np-2; ip++ ) {
    Double_t ang = dang*(ip-2)+delta;
    xp[ip] = -radhole*cos(ang);    ;
    yp[ip] = holevtx - radhole*sin(ang);
  }
  
  TGeoXtru* mdlShape = new TGeoXtru(2);
  mdlShape -> SetName(shpName.c_str());
  
  mdlShape-> DefinePolygon(np,xp,yp);
  
  mdlShape -> DefineSection(0,-halfthick);
  mdlShape -> DefineSection(1,+halfthick);
  
  return mdlShape;
  
} 

//_____________________________________________________________________________
TGeoShape* GeoScintMdlVolume::BuildNearLmXtru(const UgliDbiTables& ugliTables,
                                              std::string shpName,
                                              bool isInner) {
 // Box with coil hole cut out from top (high tpos) side, with hole
 // centered on lpos = 0

  UpdateGlobalManager();

  Float_t m01 =0; // slope on lo lpos side
  Float_t m23 = 0; // slope on hi lpos side
  Build4Corners(ugliTables,m01,m23); // stores 4 corners of trap in fXY

  const int np = 15; // must be greater than 5
  Double_t xp[np] = {0}, yp[np] = {0};

  Double_t scale = Geo::GetScale(fGeoGeometry->GetAppType());
  // 26 cm radius to circle square coil hole
  Double_t halfthick = 0.5*fThick;
  Double_t skinthick = 0;
  Double_t radhole = Geo::kNearModuleHoleRad*scale;
  Double_t delta = 0;
  Double_t dang = TMath::Pi()/(np-5);
  radhole /= cos(0.5*dang); // so that polygon circumscribes hole
  
  if ( isInner ) {
    halfthick -= fSkinThick;
    skinthick = fSkinThick;
    radhole += fSkinThick;
    delta = TMath::ASin(skinthick/radhole); // small angle approx
    dang = (TMath::Pi()-2.*delta)/(np-5);
  }
    
  xp[0] = fXY[0][0] + skinthick; yp[0] = fXY[0][1] + skinthick;
  xp[1] = fXY[1][0] + skinthick; yp[1] = fXY[1][1] - skinthick;
  xp[np-2] = fXY[2][0] - skinthick; yp[np-2] = fXY[2][1] - skinthick;
  xp[np-1] = fXY[3][0] - skinthick; yp[np-1] = fXY[3][1] + skinthick;
  if ( isInner && TMath::Abs(fXY[2][0] - fXY[3][0]) > 0.0001*scale) {
    Double_t m23 = (fXY[2][1]-fXY[3][1])/(fXY[2][0]-fXY[3][0]);
    Double_t b23 = fXY[2][1] - m23*fXY[2][0];
    Double_t arctanm = atan(TMath::Abs(m23));
    Double_t delta = fSkinThick/sin(arctanm);
    xp[np-2] = (yp[np-2] - b23)/m23 - delta;
    xp[np-1] = (yp[np-1] - b23)/m23 - delta;
  }

  Double_t holevtx = fXY[2][1];
  for ( int ip = 2; ip < np-2; ip++ ) {
    Double_t ang = dang*(ip-2)+delta;
    xp[ip] = -radhole*cos(ang);    ;
    yp[ip] = holevtx - radhole*sin(ang);
  }
  
  TGeoXtru* mdlShape = new TGeoXtru(2);
  mdlShape -> SetName(shpName.c_str());
  
  mdlShape-> DefinePolygon(np,xp,yp);
  
  mdlShape -> DefineSection(0,-halfthick);
  mdlShape -> DefineSection(1,+halfthick);
  
  return mdlShape;

} 

//_____________________________________________________________________________
TGeoShape* GeoScintMdlVolume::BuildNearlMXtru(const UgliDbiTables& ugliTables,
                                              std::string shpName,
                                              bool isInner) {
 // Box with coil hole cut out from bottom (lo tpos) side, with hole
 // centered on lpos = 0

  UpdateGlobalManager();

  Float_t m01 =0; // slope on lo lpos side
  Float_t m23 = 0; // slope on hi lpos side
  Build4Corners(ugliTables,m01,m23); // stores 4 corners of trap in fXY

  const int np = 15; // must be greater than 5
  Double_t xp[np] = {0}, yp[np] = {0};

  Double_t scale = Geo::GetScale(fGeoGeometry->GetAppType());
  // 26 cm radius to circle square coil hole
  Double_t halfthick = 0.5*fThick;
  Double_t skinthick = 0;
  Double_t radhole = Geo::kNearModuleHoleRad*scale;
  Double_t delta = 0;
  Double_t dang = TMath::Pi()/(np-5);
  radhole /= cos(0.5*dang); // so that polygon circumscribes hole
  
  if ( isInner ) {
    halfthick -= fSkinThick;
    skinthick = fSkinThick;
    radhole += fSkinThick;
    delta = TMath::ASin(skinthick/radhole); // small angle approx
    dang = (TMath::Pi()-2.*delta)/(np-5);
  }
    
  xp[0] = fXY[0][0] + skinthick; yp[0] = fXY[0][1] + skinthick;
  xp[1] = fXY[1][0] + skinthick; yp[1] = fXY[1][1] - skinthick;
  xp[2] = fXY[2][0] - skinthick; yp[2] = fXY[2][1] - skinthick;
  xp[3] = fXY[3][0] - skinthick; yp[3] = fXY[3][1] + skinthick;
  if ( isInner && TMath::Abs(fXY[2][0] - fXY[3][0]) > 0.0001*scale) {
    Double_t m23 = (fXY[2][1]-fXY[3][1])/(fXY[2][0]-fXY[3][0]);
    Double_t b23 = fXY[2][1] - m23*fXY[2][0];
    Double_t arctanm = atan(TMath::Abs(m23));
    Double_t delta = fSkinThick/sin(arctanm);
    xp[2] = (yp[2] - b23)/m23 - delta;
    xp[3] = (yp[3] - b23)/m23 - delta;
  }

  Double_t holevtx = fXY[3][1];
  for ( int ip = 4; ip < np; ip++ ) {
    Double_t ang = dang*(ip-4)+delta;
    xp[ip] = radhole*cos(ang);    ;
    yp[ip] = holevtx + radhole*sin(ang);
  }
  
  TGeoXtru* mdlShape = new TGeoXtru(2);
  mdlShape -> SetName(shpName.c_str());
  
  mdlShape-> DefinePolygon(np,xp,yp);
  
  mdlShape -> DefineSection(0,-halfthick);
  mdlShape -> DefineSection(1,+halfthick);
  
  return mdlShape;

} 

//_____________________________________________________________________________
TGeoShape* GeoScintMdlVolume::BuildNearIXtru(const UgliDbiTables& ugliTables,
                                             std::string shpName,
                                             bool isInner) {
 // Box with coil hole cut out from lo lpos, hi tpos corner

  UpdateGlobalManager();
  
  Float_t m01 =0; // slope on lo lpos side
  Float_t m23 = 0; // slope on hi lpos side
  Build4Corners(ugliTables,m01,m23); // stores 4 corners of trap in fXY

  // For module type 'I', strips 48-56 are long, strips 57-67 are short
  // To define corner cutout, need tpos of hi tpos side of last long strip(56)
  PlexStripEndId seid(fScintMdlId,56);
  const UgliDbiStrip* scStrip = ugliTables.GetDbiStripById(seid);
  std::string stpName = GeoGeometry::GetGeoCompatibleName(seid.AsString("s"));
  const GeoStripVolume* stpVol = dynamic_cast<const GeoStripVolume*>
                            (gGeoManager->GetVolume(stpName.c_str()));
  if ( !stpVol ) {
    MSG("Geo",Msg::kFatal) 
       << "Unable to find pre-built stripvol for volume "
       << stpName.c_str() << endl;
    abort();
  }
  Float_t ypositype = GetScale()*(scStrip->GetTPosRelMdl())
                     + stpVol->GetHalfWidth() + fSkinThick;
  
  const int np = 6; 
  Double_t xp[np] = {0}, yp[np] = {0};

  Double_t halfthick = 0.5*fThick;
  Double_t skinthick = 0;
  
  if ( isInner ) {
    halfthick -= fSkinThick;
    skinthick = fSkinThick;
  }
    
  xp[0] = fXY[0][0] + skinthick; yp[0] = fXY[0][1] + skinthick;
  xp[1] = fXY[0][0] + skinthick; yp[1] = ypositype - skinthick;
  xp[2] = fXY[1][0] + skinthick; yp[2] = ypositype - skinthick;
  xp[3] = fXY[1][0] + skinthick; yp[3] = fXY[1][1] - skinthick;
  xp[4] = fXY[2][0] - skinthick; yp[4] = fXY[2][1] - skinthick;
  xp[5] = fXY[3][0] - skinthick; yp[5] = fXY[3][1] + skinthick;
   
  TGeoXtru* mdlShape = new TGeoXtru(2);
  mdlShape -> SetName(shpName.c_str());
  
  mdlShape-> DefinePolygon(np,xp,yp);
  
  mdlShape -> DefineSection(0,-halfthick);
  mdlShape -> DefineSection(1,+halfthick);
  
  return mdlShape;

} 

//_____________________________________________________________________________
TGeoShape* GeoScintMdlVolume::BuildNeariXtru(const UgliDbiTables& ugliTables,
                                             std::string shpName,
                                             bool isInner) {
 // Box with coil hole cut out from lo lpos, lo tpos corner

  UpdateGlobalManager();

  Float_t m01 =0; // slope on lo lpos side
  Float_t m23 = 0; // slope on hi lpos side
  Build4Corners(ugliTables,m01,m23); // stores 4 corners of trap in fXY

  // For module type 'i', strips 0-10 are short, strips 11-19 are long
  // To define corner cutout, need tpos of lo tpos side of first long strip(11)
  PlexStripEndId seid(fScintMdlId,11);
  const UgliDbiStrip* scStrip = ugliTables.GetDbiStripById(seid);
  std::string stpName = GeoGeometry::GetGeoCompatibleName(seid.AsString("s"));
  const GeoStripVolume* stpVol = dynamic_cast<const GeoStripVolume*>
                            (gGeoManager->GetVolume(stpName.c_str()));

  if ( !stpVol ) {
    MSG("Geo",Msg::kFatal) 
       << "Unable to find pre-built stripvol for volume "
       << stpName.c_str() << endl;
    abort();
  }
  Float_t ypositype = GetScale()*(scStrip->GetTPosRelMdl())
                    - stpVol->GetHalfWidth() - fSkinThick;

  const int np = 6; 
  Double_t xp[np] = {0}, yp[np] = {0};

  Double_t halfthick = 0.5*fThick;
  Double_t skinthick = 0;
  
  if ( isInner ) {
    halfthick -= fSkinThick;
    skinthick = fSkinThick;
  }
    
  xp[0] = fXY[0][0] + skinthick; yp[0] = fXY[0][1] + skinthick;
  xp[1] = fXY[0][0] + skinthick; yp[1] = ypositype + skinthick;
  xp[2] = fXY[1][0] + skinthick; yp[2] = ypositype + skinthick;
  xp[3] = fXY[1][0] + skinthick; yp[3] = fXY[1][1] - skinthick;
  xp[4] = fXY[2][0] - skinthick; yp[4] = fXY[2][1] - skinthick;
  xp[5] = fXY[3][0] - skinthick; yp[5] = fXY[3][1] + skinthick;
   
  TGeoXtru* mdlShape = new TGeoXtru(2);
  mdlShape -> SetName(shpName.c_str());
  
  mdlShape-> DefinePolygon(np,xp,yp);
  
  mdlShape -> DefineSection(0,-halfthick);
  mdlShape -> DefineSection(1,+halfthick);
  
  return mdlShape;

} 

//_____________________________________________________________________________
TGeoShape* GeoScintMdlVolume::BuildNearXtru(const UgliDbiTables& ugliTables,
                                            std::string shpName, bool isInner){
 // Box with coil hole cut out from bottom (lo tpos) side, with hole
 // centered on lpos = 0

  UpdateGlobalManager();

  Float_t m01 =0; // slope on lo lpos side
  Float_t m23 = 0; // slope on hi lpos side
  Build4Corners(ugliTables,m01,m23); // stores 4 corners of trap in fXY

  const int np = 4; 
  Double_t xp[np] = {0}, yp[np] = {0};

  Double_t halfthick = 0.5*fThick;
  Double_t skinthick = 0;
  
  if ( isInner ) {
    halfthick -= fSkinThick;
    skinthick = fSkinThick;
  }
    
  xp[0] = fXY[0][0] + skinthick; yp[0] = fXY[0][1] + skinthick;
  xp[1] = fXY[1][0] + skinthick; yp[1] = fXY[1][1] - skinthick;
  xp[2] = fXY[2][0] - skinthick; yp[2] = fXY[2][1] - skinthick;
  xp[3] = fXY[3][0] - skinthick; yp[3] = fXY[3][1] + skinthick;

  if ( isInner && m23 > 0.0001) {
    Double_t b23 = fXY[2][1] - m23*fXY[2][0];
    Double_t arctanm = atan(TMath::Abs(m23));
    Double_t delta = fSkinThick/sin(arctanm);
    xp[2] = (yp[2] - b23)/m23 - delta;
    xp[3] = (yp[3] - b23)/m23 - delta;
  }

  if ( isInner && m01 > 0.0001) {
    Double_t b01 = fXY[1][1] - m01*fXY[1][0];
    Double_t arctanm = atan(TMath::Abs(m01));
    Double_t delta = fSkinThick/sin(arctanm);
    xp[1] = (yp[1] - b01)/m01 + delta;
    xp[0] = (yp[0] - b01)/m01 + delta;
  }
  
  TGeoXtru* mdlShape = new TGeoXtru(2);
  mdlShape -> SetName(shpName.c_str());
  
  mdlShape-> DefinePolygon(np,xp,yp);
  
  mdlShape -> DefineSection(0,-halfthick);
  mdlShape -> DefineSection(1,+halfthick);
  
  return mdlShape;

} 

//_____________________________________________________________________________
TGeoShape* GeoScintMdlVolume::BuildFarDFXtru(const UgliDbiTables& ugliTables,
                                             std::string shpName,
                                             bool isInner) {
 // Box with coil hole cut out from bottom

  UpdateGlobalManager();

  Float_t m01 =0; // slope on lo lpos side
  Float_t m23 = 0; // slope on hi lpos side
  Build4Corners(ugliTables,m01,m23); // stores 4 corners of trap in fXY

  const int np = 15; // greater than 5
  Double_t xp[np] = {0}, yp[np] = {0};

  Double_t halfthick = 0.5*fThick;
  Double_t skinthick = 0;
  Double_t radhole = Geo::kFarModuleHoleRad*GetScale();
  Double_t delta = 0;
  Double_t dang = TMath::Pi()/(np-5);
  radhole /= cos(0.5*dang);
  
  if ( isInner ) {
    halfthick -= fSkinThick;
    skinthick = fSkinThick;
    radhole += fSkinThick;
    delta = TMath::ASin(skinthick/radhole); // small angle approx
    dang = (TMath::Pi()-2.*delta)/(np-5);
  }
    
  xp[0] = fXY[0][0] + skinthick; yp[0] = fXY[0][1] + skinthick;
  xp[1] = fXY[1][0] + skinthick; yp[1] = fXY[1][1] - skinthick;
  xp[2] = fXY[2][0] - skinthick; yp[2] = fXY[2][1] - skinthick;
  xp[3] = fXY[3][0] - skinthick; yp[3] = fXY[3][1] + skinthick;
  Double_t holevtx = fXY[3][1];
  for ( int ip = 4; ip < np; ip++ ) {
    Double_t ang = dang*(ip-4)+delta;
    xp[ip] = radhole*cos(ang);    ;
    yp[ip] = holevtx + radhole*sin(ang);
  }
  
  TGeoXtru* mdlShape = new TGeoXtru(2);
  mdlShape -> SetName(shpName.c_str());
  
  mdlShape-> DefinePolygon(np,xp,yp);
  
  mdlShape -> DefineSection(0,-halfthick);
  mdlShape -> DefineSection(1,+halfthick);
  
  return mdlShape;

} 

//_____________________________________________________________________________
void GeoScintMdlVolume::Build4Corners(const UgliDbiTables& ugliTables,
                                      Float_t& m01, Float_t& m23) {
 // Private routine to help with build of module shape
 // Fills array fXY[4][2] with xy positions defining the 4 corners of 
 // the trapezoid, and returns (through reference arguments) the slope of the
 // lo lpos (m01) and hi lpos (m23) sides.  The ordering of the 4 corners
 // is clockwise from the lo lpos, lo tpos position, e.g. 0)(lolpos,lotpos),
 // 1)(lolpos,hitpos), 2)(hilpos,hitpos), 3)(hilpos,lotpos)  
 // Assumes that strip volumes have already been built.

  UpdateGlobalManager();

  m01 = 0; m23 = 0;  // null slope implies vertical
  Float_t scale = GetScale();
  
  UgliDbiStructHash modulestructhash(fScintMdlId);
  const UgliDbiScintMdlStruct* scModuleStruct 
       = ugliTables.fScintMdlStructTbl
        .GetRowByIndex(modulestructhash.HashAsScintMdl());

  // Lo and hi strips associated with this module
  Int_t lostrip = scModuleStruct -> GetFirstStrip();
  Int_t histrip = scModuleStruct -> GetLastStrip(); 
  
  // The assumption made to determine the boundaries of a given module
  // is that zrot of strips relative to module is always null, and that
  // lpos of strips relative to module follow the boundaries of the trapezoid
  // defined by the 4 lo/hi strip corners. Alignment variations in tpos
  // are allowed.

  // Need lo strip position to set lower corner boundaries
  PlexStripEndId loseid(fScintMdlId,lostrip);
  const UgliDbiStrip* scStrip = ugliTables.GetDbiStripById(loseid);
  std::string stpName=GeoGeometry::GetGeoCompatibleName(loseid.AsString("s"));
  const GeoStripVolume* stpVol = dynamic_cast<const GeoStripVolume*>
              (gGeoManager -> GetVolume(stpName.c_str()));
  if ( !stpVol ) {
    MSG("Geo",Msg::kFatal) 
       << "Unable to find pre-built stripvol for volume "
       << stpName.c_str() << endl;
    abort();
  }
  
  Float_t lostriphalfwidth = stpVol->GetHalfWidth() + fSkinThick;
  
  // lo lpos, hi tpos corner of lo strip
  fXY[0][0] = scale*(scStrip->GetLPosRelMdl())-stpVol->GetHalfLength()
             -fSkinThick;
  fXY[0][1] = scale*(scStrip->GetTPosRelMdl())+stpVol->GetHalfWidth()
             +fSkinThick;
    
  // hi lpos, hi tpos corner of lo strip
  fXY[3][0] = scale*(scStrip->GetLPosRelMdl())+stpVol->GetHalfLength()
             +fSkinThick;
  fXY[3][1] = scale*(scStrip->GetTPosRelMdl())+stpVol->GetHalfWidth()
             +fSkinThick;

  // Need hi strip tpos first to set upper corner boundaries
  PlexStripEndId hiseid(fScintMdlId,histrip);
  scStrip = ugliTables.GetDbiStripById(hiseid);
  stpName=GeoGeometry::GetGeoCompatibleName(hiseid.AsString("s"));

  stpVol = dynamic_cast<GeoStripVolume*>
           (gGeoManager -> GetVolume(stpName.c_str()));
  if ( !stpVol ) {
    MSG("Geo",Msg::kFatal) 
       << "Unable to find pre-built stripvol for volume "
       << stpName.c_str() << endl;
    abort();
  }

  Double_t histriphalfwidth = stpVol->GetHalfWidth() + fSkinThick;

  // lo lpos, hi tpos of histrip 
  fXY[1][0] = scale*(scStrip->GetLPosRelMdl())-stpVol->GetHalfLength()
             -fSkinThick;
  fXY[1][1] = scale*(scStrip->GetTPosRelMdl())+stpVol->GetHalfWidth()
             +fSkinThick;

  // hi lpos, hi tpos of hi strip
  fXY[2][0] = scale*(scStrip->GetLPosRelMdl())+stpVol->GetHalfLength()
             +fSkinThick;
  fXY[2][1] = scale*(scStrip->GetTPosRelMdl())+stpVol->GetHalfWidth()
             +fSkinThick;

  // Calculate slope on p0/p1 side
  UChar_t mdlType = fScintMdlId.GetModuleType();
  if ( mdlType != 'i' && mdlType != 'I' ) {
    // Not for near detector modules with corner cut-outs. Fix me.
    if ( TMath::Abs(fXY[1][0] - fXY[0][0]) >= 0.0001*scale ) {
      m01 = (fXY[1][1] - fXY[0][1])/
            (fXY[1][0] - fXY[0][0]);
      if ( m01 > 0 ) {
        fXY[0][0] -= 2.*lostriphalfwidth/m01; // pos slope
      }
      else {
        fXY[1][0] += 2.*histriphalfwidth/m01; // neg slope
      }
    }
  }
  
  // Calculate slope on p3/p2 side
  if ( TMath::Abs(fXY[2][0] - fXY[3][0]) >= 0.0001*scale ) {
    m23 = (fXY[2][1] - fXY[3][1])/
          (fXY[2][0] - fXY[3][0]);
    if ( m23 < 0 ) {
      fXY[3][0] -= 2.*lostriphalfwidth/m23; // neg slope
    }
    else {
      fXY[2][0] += 2.*histriphalfwidth/m23; // pos slope
    }
  }

  fXY[0][1] -= 2.*lostriphalfwidth;
  fXY[3][1] -= 2.*lostriphalfwidth;

  //for (int ic = 0; ic < 4; ic++ ) {
  //  cout << ic << ")" << fXY[ic][0] << "," << fXY[ic][1] << endl;
  //}
  
}

//_____________________________________________________________________________
void GeoScintMdlVolume::BuildHalfY(const UgliDbiTables& ugliTables) {
 // Private routine to help with build of module shape
 // Fills fXY[0][1] and fXY[1][1] matrix elements only to be used
 // in building of scaled shapes.

  UpdateGlobalManager();

  Float_t scale = GetScale();

  UgliDbiStructHash modulestructhash(fScintMdlId);
  const UgliDbiScintMdlStruct* scModuleStruct
       = ugliTables.fScintMdlStructTbl
        .GetRowByIndex(modulestructhash.HashAsScintMdl());

  // Lo and hi strips associated with this module
  Int_t lostrip = scModuleStruct -> GetFirstStrip();
  Int_t histrip = scModuleStruct -> GetLastStrip();

  // Need lo strip position to set lower corner boundaries
  PlexStripEndId loseid(fScintMdlId,lostrip);
  const UgliDbiStrip* scStrip = ugliTables.GetDbiStripById(loseid);

  Float_t striphalfwidth = 0.5*(Geo::kStripWidth)*scale + fSkinThick;

  fXY[0][1] = scale*(scStrip->GetTPosRelMdl())-striphalfwidth;

  // Need hi strip tpos to set upper corner boundaries
  PlexStripEndId hiseid(fScintMdlId,histrip);
  scStrip = ugliTables.GetDbiStripById(hiseid);

  // lo lpos, hi tpos of histrip
  fXY[1][1] = scale*(scStrip->GetTPosRelMdl())+striphalfwidth;

}

//_____________________________________________________________________________
Int_t GeoScintMdlVolume::NumberOfStrips() const {
 // Return number of strips in volume
  
  UpdateGlobalManager();

  return (this -> GetAirNode() -> GetVolume() -> GetNdaughters());
  
}

---