Geometry.cxx

Go to the documentation of this file.

#include "Geometry/Geometry.h"
extern "C" {
#include <sys/types.h>
#include <sys/stat.h>
}
#include <cassert>
#include <map>
#include <iostream>
#include "TObjArray.h"
#include "TGeoManager.h"
#include "TGeoNode.h"
#include "TGeoBBox.h"
#include "TMath.h"
#include "Geometry/Exception.h"
#include "Geometry/PlaneGeo.h"
#include "Geometry/CellGeo.h"
#include "Geometry/CellUniqueId.h"
#include "JobControl/Exception.h"
#include "RawData/DAQHeader.h"
using namespace geo;

Geometry* Geometry::fInstance = 0;

static bool plane_sort(const PlaneGeo* p1, const PlaneGeo* p2) 
{
  double xyz1[3], xyz2[3];
  p1->Cell(0).GetCenter(xyz1);
  p2->Cell(0).GetCenter(xyz2);
  return xyz1[2]<xyz2[2];
}

//......................................................................

Geometry& Geometry::Instance(int runnum, short detId)
{
  int oldrunnum = -1;
  int olddetid  = -1;
  // If our current geometry doesn't match request, get rid of it
  if (fInstance!=0) {
    if (fInstance->fRunNumber != runnum || fInstance->fDetId != detId) {
      oldrunnum = fInstance->fRunNumber;
      olddetid  = fInstance->fDetId;
      delete fInstance;
      fInstance = 0;
      gGeoManager->Clear();
    }
  }

  // Build a new geometry if we need to
  if (fInstance==0) {
    static char* fname[3] = {"neardet.gdml","fardet.gdml","ipnd.gdml"};
    if (detId!=rawdata::kNear &&
        detId!=rawdata::kFar  &&
        detId!=rawdata::kIPND) {
      char msg[256];
      sprintf(msg,"Unknown detector id number: %d",(int)detId);
      throw Exception(__FILE__,__LINE__,Exception::BAD_GEO_CONFIG,msg);
    }
    if (runnum<0) {
      char msg[256];
      std::sprintf(msg,"Bad run number: %d",runnum);
      throw Exception(__FILE__,__LINE__,Exception::BAD_GEO_CONFIG,msg);
    }
    fInstance = new Geometry(fname[detId]);
    std::cout << "Geometry updated from "
              << olddetid << ":" << oldrunnum << " to "
              << detId    << ":" << runnum    << std::endl;
    fInstance->fDetId     = detId;
    fInstance->fRunNumber = runnum;
  }
  return *fInstance;
}

//......................................................................

Geometry::Geometry(const char* gdmlfile) 
{
  std::string fname;

  // Test if the gdml file can be openned using exactly the name we've
  // been passed
  while (1) {
    struct stat sb;
    fname = gdmlfile;
    if (stat(fname.c_str(), &sb)==0) break;

    const char* srt_private = getenv("SRT_PRIVATE_CONTEXT");
    if (srt_private!=0) {
      fname = srt_private;
      fname += "/Geometry/gdml/";
      fname += gdmlfile;
      if (stat(fname.c_str(), &sb)==0) break;
    }

    const char* srt_public = getenv("SRT_PUBLIC_CONTEXT");
    if (srt_public!=0) {
      fname = srt_public;
      fname += "/Geometry/gdml/";
      fname += gdmlfile;
      if (stat(fname.c_str(), &sb)==0) break;
    }
    
    // Failed to resolve the file name
    throw Exception(__FILE__, __LINE__, Exception::FILE_NOT_FOUND,
                    "GDML file not found.");
  }

  // Open the GDML file
  TGeoManager::Import(fname.c_str());
  this->SetDrawOptions();
  
  std::vector<const TGeoNode*> n(16);
  n[0] = gGeoManager->GetTopNode();
  this->FindPlanes(n, 0);
  sort(fPlanes.begin(), fPlanes.end(), plane_sort);
  
  this->SetDetectorSize();
  this->BuildMaps();
}

//......................................................................
TGeoManager* Geometry::ROOTGeoManager() const
{
  return gGeoManager;
}

//......................................................................

Geometry::~Geometry() 
{
  for (unsigned int i=0; i<fPlanes.size(); ++i) {
    if (fPlanes[i]) { delete fPlanes[i]; fPlanes[i] = 0; }
  }
}

//......................................................................
int Geometry::CurrentCell(int* planeid, int* cellid) const 
{
  try {
    this->IdToCell(this->CurrentCellId(), planeid, cellid);
  }
  catch (Exception e) {
    *planeid = -1;
    *cellid  = -1;
    return -1;
  }
  return 1;
}

//......................................................................
const CellUniqueId Geometry::CurrentCellId() const
{
  return PathToUniqueId(gGeoManager->GetPath());
}

//......................................................................

const TGeoMaterial* Geometry::Material(double x, double y, double z) const 
{
  const TGeoNode* node = gGeoManager->FindNode(x,y,z);
  return node->GetMedium()->GetMaterial();
}

//......................................................................

double Geometry::DetHalfWidth()  const 
{
  return fDetHalfWidth;
}

//......................................................................

double Geometry::DetHalfHeight() const 
{
  return fDetHalfHeight;
}

//......................................................................

double Geometry::DetLength() const
{ 
  return fDetLength;
}

//......................................................................

double Geometry::SurfaceY() const 
{
  switch (fDetId) {
  case 0: return 130.0E2; // near det: surface ~130m up from det. center
  case 1: return 9.94E2;  // far det: surface ~10m up from det. center
  case 2: return -1.966E2;  // IPND: surface ~2 meter below det, center
  default: throw Exception(__FILE__,__LINE__,Exception::BAD_GEO_CONFIG);
  }
  return 0.0;
}

void Geometry::WorldBox(double* xlo, double* xhi,
                        double* ylo, double* yhi,
                        double* zlo, double* zhi) const
{
  const TGeoShape* s = gGeoManager->GetVolume("vWorld")->GetShape();
  assert(s);

  if (xlo || xhi) {
    double x1, x2;
    s->GetAxisRange(1,x1,x2); if (xlo) *xlo = x1; if (xhi) *xhi = x2;
  }
  if (ylo || yhi) {
    double y1, y2;
    s->GetAxisRange(2,y1,y2); if (ylo) *ylo = y1; if (yhi) *yhi = y2;
  }
  if (zlo || zhi) {
    double z1, z2;
    s->GetAxisRange(3,z1,z2); if (zlo) *zlo = z1; if (zhi) *zhi = z2;
  }
}

//......................................................................

void Geometry::SetDrawOptions() { }

//......................................................................

void Geometry::FindPlanes(std::vector<const TGeoNode*>& n, 
                          unsigned int depth)
{
  const char* nm = n[depth]->GetName();
  if (strncmp(nm,"vPlane", 6)==0) {
    this->MakePlane(n, depth);
    return;
  }

  // Explore the next layer down
  unsigned int deeper = depth+1;
  if (deeper>=n.size()) {
    throw Exception(__FILE__,__LINE__,Exception::BAD_NODE,
                    "Exceeded maximum depth.");
  }
  const TGeoVolume* v = n[depth]->GetVolume();
  int nd = v->GetNdaughters();
  for (int i=0; i<nd; ++i) {
    n[deeper] = v->GetNode(i);
    this->FindPlanes(n, deeper);
  }
}

//......................................................................

void Geometry::MakePlane(std::vector<const TGeoNode*>& n,
                         unsigned int depth)
{
  fPlanes.push_back(new PlaneGeo(n, depth));
}

//......................................................................
const PlaneGeo& Geometry::Plane(unsigned int i) const 
{
  if (i>=fPlanes.size()) throw Exception(__FILE__,__LINE__,
                                         Exception::BAD_PLANE_NUMBER,
                                         "Plane index out of range");
  return *fPlanes[i];
}

//......................................................................
const CellGeo& Geometry::IdToCell(const CellUniqueId& id, 
                                  int* iplane, 
                                  int* icell) const 
{
  UIDMap::const_iterator itr = fIdMap.find(id);
  if (itr == fIdMap.end()) {
    throw Exception(__FILE__,__LINE__,Exception::BAD_CELL_UNIQUE_ID);
  }
  *iplane = itr->second.first;
  *icell  = itr->second.second;
  return this->Plane(*iplane).Cell(*icell);
}

//......................................................................
const std::set<unsigned int>& Geometry::GetPlanesByView(View_t v) const 
{
  // Return the user's choice, fall back on all planes
  if      (v==kX) return fXplanes;
  else if (v==kY) return fYplanes;
  return fAllPlanes;
}

//......................................................................
const unsigned int Geometry::NextPlaneInView(unsigned int p1, int d) const 
{
  assert(p1<fPlanes.size());

  // March along from p1 in the direction d until we find a match
  int s = (d>0 ? 1 : -1);
  if (s<0 && p1==0)                return kPLANE_NOT_FOUND;
  if (s>0 && p1==fPlanes.size()-1) return kPLANE_NOT_FOUND;
  for (unsigned int p2=p1+d; 1; p2+=s) {
    if (fPlanes[p1]->View() == fPlanes[p2]->View()) return p2;
    if (p2==fPlanes.size()-1) return kPLANE_NOT_FOUND;
    if (p2==0)                return kPLANE_NOT_FOUND;
  }
}

//......................................................................
const unsigned int Geometry::NextPlaneOtherView(unsigned int p1, int d) const 
{
  assert(p1<fPlanes.size());

  // March along from p1 in the direction d until we find a match
  int s = (d>0 ? 1 : -1);
  if (s<0 && p1==0)                return kPLANE_NOT_FOUND;
  if (s>0 && p1==fPlanes.size()-1) return kPLANE_NOT_FOUND;
  for (unsigned int p2=p1+d; 1; p2+=s) {
    if (fPlanes[p1]->View() != fPlanes[p2]->View()) return p2;
    if (p2==fPlanes.size()) return kPLANE_NOT_FOUND;
    if (p2==0)              return kPLANE_NOT_FOUND;
  }
}

//......................................................................
void Geometry::BuildMaps ()
{
  // Sets which list planes by view
  for (unsigned int i=0; i<fPlanes.size(); ++i) {
    fAllPlanes.insert(i);
    if      (fPlanes[i]->View()==kX) fXplanes.insert(i);
    else if (fPlanes[i]->View()==kY) fYplanes.insert(i);
    else {
      throw Exception(__FILE__,__LINE__,Exception::BAD_GEO_CONFIG,
                      "Bad plane view!");
    }
  }
  
  // Unique Id -> cell map
  for (unsigned int i=0; i<fPlanes.size(); ++i) {
    const PlaneGeo* p = fPlanes[i];
    for (int j=0; j<p->Ncells(); ++j) {
      // Assert that the cell ids are in fact unique
      if (fIdMap.find(this->Plane(i).Cell(j).Id())!=fIdMap.end()) {
        throw Exception(__FILE__,__LINE__,Exception::BAD_CELL_UNIQUE_ID,
                        "Duplicate Cell UIDs");
      }
      PlaneCellPair p(i,j);
      fIdMap[this->Plane(i).Cell(j).Id()] = p;
    }
  }
}

//......................................................................
double Geometry::TotalMass() const
{
  double mass = 0.;

  const TGeoVolume *enclosure = gGeoManager->FindVolumeFast("vDetEnclosure");
  
  for(int d = 0; d < enclosure->GetNdaughters(); ++d){
    const TGeoNode* node = enclosure->GetNode(d);
    const char* nm = node->GetName();
    if(strncmp(nm,"vSuperBlock",11) == 0){
      std::string name = nm;
      name.erase(name.find('_'));
      mass += gGeoManager->FindVolumeFast(name.c_str())->Weight();
    }
  }
  return mass;
}

//......................................................................
double Geometry::MassBetweenPoints(double *p1, double *p2) const
{

  double columnD = 0.;

  double length = TMath::Sqrt(TMath::Power(p2[0]-p1[0], 2.)
                              + TMath::Power(p2[1]-p1[1], 2.)
                              + TMath::Power(p2[2]-p1[2], 2.));
  double dxyz[3] = {(p2[0]-p1[0])/length, (p2[1]-p1[1])/length, (p2[2]-p1[2])/length}; 

  gGeoManager->InitTrack(p1,dxyz);

  TGeoNode *node = gGeoManager->GetCurrentNode();

  while(!gGeoManager->IsSameLocation(p2[0], p2[1], p2[2])){
    gGeoManager->FindNextBoundary();
    columnD += gGeoManager->GetStep()*node->GetMedium()->GetMaterial()->GetDensity();
    
    node = gGeoManager->Step();
  }//end loop to get to volume of second point

  const double *current = gGeoManager->GetCurrentPoint();
  length = TMath::Sqrt(TMath::Power(p2[0]-current[0], 2.)
                       + TMath::Power(p2[1]-current[1], 2.)
                       + TMath::Power(p2[2]-current[2], 2.));
  columnD += length*node->GetMedium()->GetMaterial()->GetDensity();

  return columnD;
}

//......................................................................

void Geometry::SetDetectorSize() 
{
  double dummy;

  // Compute the height and width based on the sizes of the vertical
  // and horizontal planes. Remember, in the local plane frame z goes
  // along the length of the modules. Hence the use of axis=3 to get
  // the detector height (from the vertical planes) and width (from
  // the horizontal planes).
  const TGeoVolume* vvp = gGeoManager->GetVolume("vPlaneV");
  const TGeoVolume* hvp = gGeoManager->GetVolume("vPlaneH");

  // Handle the case of the IPND where the volumes are names slightly
  // different accounting for the different widths of the modules.
  if (vvp==0) vvp = gGeoManager->GetVolume("vPlaneVIPND");
  if (hvp==0) hvp = gGeoManager->GetVolume("vPlaneHIPND");
  if (vvp==0 && hvp==0) {
    vvp = gGeoManager->GetVolume("vPlaneVIPND");
    hvp = gGeoManager->GetVolume("vPlaneHIPND");
  }
  if (vvp==0 || hvp==0) {
    throw Exception(__FILE__,__LINE__,
                    Exception::BAD_GEO_CONFIG,
                    "Unable to find shapes to set detector size");
  }
  
  const TGeoShape* vp = vvp->GetShape();
  const TGeoShape* hp = hvp->GetShape();
  vp->GetAxisRange(3,dummy,fDetHalfHeight);
  hp->GetAxisRange(3,dummy,fDetHalfWidth);

  // Remember -- in the local plane coordinate frame, the "depth" is
  // along y, hence axis=2. In the world coordinate this is along z.
  double vplanez1, vplanez2;
  vp->GetAxisRange(2,vplanez1,vplanez2);
  
  // Compute the detector length based on the center positions of
  // cells in the first and last plane. Adjust to account for the
  // plane widths. This ajustment works if the vertical and horizontal
  // planes have the same depth, or if the detector begins and ends
  // with vertical planes. As of the TDR both of these conditions was
  // true for all the NOvA detectors.
  double xyz1[3] = {0,0,0};
  double xyz2[3] = {0,0,0};
  this->Plane(0).               Cell(0).GetCenter(xyz1);
  this->Plane(fPlanes.size()-1).Cell(0).GetCenter(xyz2);
  fDetLength = (xyz2[2]-xyz1[2])+(vplanez2-vplanez1);
}

---