00001
00002
00003
00004
00005
00006
00007 #include "Geometry/Geometry.h"
00008 extern "C" {
00009 #include <sys/types.h>
00010 #include <sys/stat.h>
00011 }
00012 #include <cassert>
00013 #include <map>
00014 #include <iostream>
00015 #include "TObjArray.h"
00016 #include "TGeoManager.h"
00017 #include "TGeoNode.h"
00018 #include "TGeoBBox.h"
00019 #include "TMath.h"
00020 #include "Geometry/Exception.h"
00021 #include "Geometry/PlaneGeo.h"
00022 #include "Geometry/CellGeo.h"
00023 #include "Geometry/CellUniqueId.h"
00024 #include "JobControl/Exception.h"
00025 #include "RawData/DAQHeader.h"
00026 using namespace geo;
00027
00028 Geometry* Geometry::fInstance = 0;
00029
00030 static bool plane_sort(const PlaneGeo* p1, const PlaneGeo* p2)
00031 {
00032 double xyz1[3], xyz2[3];
00033 p1->Cell(0).GetCenter(xyz1);
00034 p2->Cell(0).GetCenter(xyz2);
00035 return xyz1[2]<xyz2[2];
00036 }
00037
00038
00039
00049 Geometry& Geometry::Instance(int runnum, short detId)
00050 {
00051 int oldrunnum = -1;
00052 int olddetid = -1;
00053
00054 if (fInstance!=0) {
00055 if (fInstance->fRunNumber != runnum || fInstance->fDetId != detId) {
00056 oldrunnum = fInstance->fRunNumber;
00057 olddetid = fInstance->fDetId;
00058 delete fInstance;
00059 fInstance = 0;
00060 gGeoManager->Clear();
00061 }
00062 }
00063
00064
00065 if (fInstance==0) {
00066 static char* fname[3] = {"neardet.gdml","fardet.gdml","ipnd.gdml"};
00067 if (detId!=rawdata::kNear &&
00068 detId!=rawdata::kFar &&
00069 detId!=rawdata::kIPND) {
00070 char msg[256];
00071 sprintf(msg,"Unknown detector id number: %d",(int)detId);
00072 throw Exception(__FILE__,__LINE__,Exception::BAD_GEO_CONFIG,msg);
00073 }
00074 if (runnum<0) {
00075 char msg[256];
00076 std::sprintf(msg,"Bad run number: %d",runnum);
00077 throw Exception(__FILE__,__LINE__,Exception::BAD_GEO_CONFIG,msg);
00078 }
00079 fInstance = new Geometry(fname[detId]);
00080 std::cout << "Geometry updated from "
00081 << olddetid << ":" << oldrunnum << " to "
00082 << detId << ":" << runnum << std::endl;
00083 fInstance->fDetId = detId;
00084 fInstance->fRunNumber = runnum;
00085 }
00086 return *fInstance;
00087 }
00088
00089
00090
00091 Geometry::Geometry(const char* gdmlfile)
00092 {
00093 std::string fname;
00094
00095
00096
00097 while (1) {
00098 struct stat sb;
00099 fname = gdmlfile;
00100 if (stat(fname.c_str(), &sb)==0) break;
00101
00102 const char* srt_private = getenv("SRT_PRIVATE_CONTEXT");
00103 if (srt_private!=0) {
00104 fname = srt_private;
00105 fname += "/Geometry/gdml/";
00106 fname += gdmlfile;
00107 if (stat(fname.c_str(), &sb)==0) break;
00108 }
00109
00110 const char* srt_public = getenv("SRT_PUBLIC_CONTEXT");
00111 if (srt_public!=0) {
00112 fname = srt_public;
00113 fname += "/Geometry/gdml/";
00114 fname += gdmlfile;
00115 if (stat(fname.c_str(), &sb)==0) break;
00116 }
00117
00118
00119 throw Exception(__FILE__, __LINE__, Exception::FILE_NOT_FOUND,
00120 "GDML file not found.");
00121 }
00122
00123
00124 TGeoManager::Import(fname.c_str());
00125 this->SetDrawOptions();
00126
00127 std::vector<const TGeoNode*> n(16);
00128 n[0] = gGeoManager->GetTopNode();
00129 this->FindPlanes(n, 0);
00130 sort(fPlanes.begin(), fPlanes.end(), plane_sort);
00131
00132 this->SetDetectorSize();
00133 this->BuildMaps();
00134 }
00135
00136
00137 TGeoManager* Geometry::ROOTGeoManager() const
00138 {
00139 return gGeoManager;
00140 }
00141
00142
00143
00144 Geometry::~Geometry()
00145 {
00146 for (unsigned int i=0; i<fPlanes.size(); ++i) {
00147 if (fPlanes[i]) { delete fPlanes[i]; fPlanes[i] = 0; }
00148 }
00149 }
00150
00151
00160 int Geometry::CurrentCell(int* planeid, int* cellid) const
00161 {
00162 try {
00163 this->IdToCell(this->CurrentCellId(), planeid, cellid);
00164 }
00165 catch (Exception e) {
00166 *planeid = -1;
00167 *cellid = -1;
00168 return -1;
00169 }
00170 return 1;
00171 }
00172
00173
00177 const CellUniqueId Geometry::CurrentCellId() const
00178 {
00179 return PathToUniqueId(gGeoManager->GetPath());
00180 }
00181
00182
00183
00184 const TGeoMaterial* Geometry::Material(double x, double y, double z) const
00185 {
00186 const TGeoNode* node = gGeoManager->FindNode(x,y,z);
00187 return node->GetMedium()->GetMaterial();
00188 }
00189
00190
00191
00192 double Geometry::DetHalfWidth() const
00193 {
00194 return fDetHalfWidth;
00195 }
00196
00197
00198
00199 double Geometry::DetHalfHeight() const
00200 {
00201 return fDetHalfHeight;
00202 }
00203
00204
00205
00206 double Geometry::DetLength() const
00207 {
00208 return fDetLength;
00209 }
00210
00211
00212
00224 double Geometry::SurfaceY() const
00225 {
00226 switch (fDetId) {
00227 case 0: return 130.0E2;
00228 case 1: return 9.94E2;
00229 case 2: return -1.966E2;
00230 default: throw Exception(__FILE__,__LINE__,Exception::BAD_GEO_CONFIG);
00231 }
00232 return 0.0;
00233 }
00234
00247 void Geometry::WorldBox(double* xlo, double* xhi,
00248 double* ylo, double* yhi,
00249 double* zlo, double* zhi) const
00250 {
00251 const TGeoShape* s = gGeoManager->GetVolume("vWorld")->GetShape();
00252 assert(s);
00253
00254 if (xlo || xhi) {
00255 double x1, x2;
00256 s->GetAxisRange(1,x1,x2); if (xlo) *xlo = x1; if (xhi) *xhi = x2;
00257 }
00258 if (ylo || yhi) {
00259 double y1, y2;
00260 s->GetAxisRange(2,y1,y2); if (ylo) *ylo = y1; if (yhi) *yhi = y2;
00261 }
00262 if (zlo || zhi) {
00263 double z1, z2;
00264 s->GetAxisRange(3,z1,z2); if (zlo) *zlo = z1; if (zhi) *zhi = z2;
00265 }
00266 }
00267
00268
00269
00270 void Geometry::SetDrawOptions() { }
00271
00272
00273
00281 void Geometry::FindPlanes(std::vector<const TGeoNode*>& n,
00282 unsigned int depth)
00283 {
00284 const char* nm = n[depth]->GetName();
00285 if (strncmp(nm,"vPlane", 6)==0) {
00286 this->MakePlane(n, depth);
00287 return;
00288 }
00289
00290
00291 unsigned int deeper = depth+1;
00292 if (deeper>=n.size()) {
00293 throw Exception(__FILE__,__LINE__,Exception::BAD_NODE,
00294 "Exceeded maximum depth.");
00295 }
00296 const TGeoVolume* v = n[depth]->GetVolume();
00297 int nd = v->GetNdaughters();
00298 for (int i=0; i<nd; ++i) {
00299 n[deeper] = v->GetNode(i);
00300 this->FindPlanes(n, deeper);
00301 }
00302 }
00303
00304
00305
00306 void Geometry::MakePlane(std::vector<const TGeoNode*>& n,
00307 unsigned int depth)
00308 {
00309 fPlanes.push_back(new PlaneGeo(n, depth));
00310 }
00311
00312
00321 const PlaneGeo& Geometry::Plane(unsigned int i) const
00322 {
00323 if (i>=fPlanes.size()) throw Exception(__FILE__,__LINE__,
00324 Exception::BAD_PLANE_NUMBER,
00325 "Plane index out of range");
00326 return *fPlanes[i];
00327 }
00328
00329
00341 const CellGeo& Geometry::IdToCell(const CellUniqueId& id,
00342 int* iplane,
00343 int* icell) const
00344 {
00345 UIDMap::const_iterator itr = fIdMap.find(id);
00346 if (itr == fIdMap.end()) {
00347 throw Exception(__FILE__,__LINE__,Exception::BAD_CELL_UNIQUE_ID);
00348 }
00349 *iplane = itr->second.first;
00350 *icell = itr->second.second;
00351 return this->Plane(*iplane).Cell(*icell);
00352 }
00353
00354
00361 const std::set<unsigned int>& Geometry::GetPlanesByView(View_t v) const
00362 {
00363
00364 if (v==kX) return fXplanes;
00365 else if (v==kY) return fYplanes;
00366 return fAllPlanes;
00367 }
00368
00369
00380 const unsigned int Geometry::NextPlaneInView(unsigned int p1, int d) const
00381 {
00382 assert(p1<fPlanes.size());
00383
00384
00385 int s = (d>0 ? 1 : -1);
00386 if (s<0 && p1==0) return kPLANE_NOT_FOUND;
00387 if (s>0 && p1==fPlanes.size()-1) return kPLANE_NOT_FOUND;
00388 for (unsigned int p2=p1+d; 1; p2+=s) {
00389 if (fPlanes[p1]->View() == fPlanes[p2]->View()) return p2;
00390 if (p2==fPlanes.size()-1) return kPLANE_NOT_FOUND;
00391 if (p2==0) return kPLANE_NOT_FOUND;
00392 }
00393 }
00394
00395
00405 const unsigned int Geometry::NextPlaneOtherView(unsigned int p1, int d) const
00406 {
00407 assert(p1<fPlanes.size());
00408
00409
00410 int s = (d>0 ? 1 : -1);
00411 if (s<0 && p1==0) return kPLANE_NOT_FOUND;
00412 if (s>0 && p1==fPlanes.size()-1) return kPLANE_NOT_FOUND;
00413 for (unsigned int p2=p1+d; 1; p2+=s) {
00414 if (fPlanes[p1]->View() != fPlanes[p2]->View()) return p2;
00415 if (p2==fPlanes.size()) return kPLANE_NOT_FOUND;
00416 if (p2==0) return kPLANE_NOT_FOUND;
00417 }
00418 }
00419
00420
00424 void Geometry::BuildMaps ()
00425 {
00426
00427 for (unsigned int i=0; i<fPlanes.size(); ++i) {
00428 fAllPlanes.insert(i);
00429 if (fPlanes[i]->View()==kX) fXplanes.insert(i);
00430 else if (fPlanes[i]->View()==kY) fYplanes.insert(i);
00431 else {
00432 throw Exception(__FILE__,__LINE__,Exception::BAD_GEO_CONFIG,
00433 "Bad plane view!");
00434 }
00435 }
00436
00437
00438 for (unsigned int i=0; i<fPlanes.size(); ++i) {
00439 const PlaneGeo* p = fPlanes[i];
00440 for (int j=0; j<p->Ncells(); ++j) {
00441
00442 if (fIdMap.find(this->Plane(i).Cell(j).Id())!=fIdMap.end()) {
00443 throw Exception(__FILE__,__LINE__,Exception::BAD_CELL_UNIQUE_ID,
00444 "Duplicate Cell UIDs");
00445 }
00446 PlaneCellPair p(i,j);
00447 fIdMap[this->Plane(i).Cell(j).Id()] = p;
00448 }
00449 }
00450 }
00451
00452
00457 double Geometry::TotalMass() const
00458 {
00459 double mass = 0.;
00460
00462 const TGeoVolume *enclosure = gGeoManager->FindVolumeFast("vDetEnclosure");
00463
00468 for(int d = 0; d < enclosure->GetNdaughters(); ++d){
00469 const TGeoNode* node = enclosure->GetNode(d);
00470 const char* nm = node->GetName();
00471 if(strncmp(nm,"vSuperBlock",11) == 0){
00472 std::string name = nm;
00473 name.erase(name.find('_'));
00474 mass += gGeoManager->FindVolumeFast(name.c_str())->Weight();
00475 }
00476 }
00477 return mass;
00478 }
00479
00480
00487 double Geometry::MassBetweenPoints(double *p1, double *p2) const
00488 {
00489
00495 double columnD = 0.;
00496
00498 double length = TMath::Sqrt(TMath::Power(p2[0]-p1[0], 2.)
00499 + TMath::Power(p2[1]-p1[1], 2.)
00500 + TMath::Power(p2[2]-p1[2], 2.));
00501 double dxyz[3] = {(p2[0]-p1[0])/length, (p2[1]-p1[1])/length, (p2[2]-p1[2])/length};
00502
00503 gGeoManager->InitTrack(p1,dxyz);
00504
00506 TGeoNode *node = gGeoManager->GetCurrentNode();
00507
00511 while(!gGeoManager->IsSameLocation(p2[0], p2[1], p2[2])){
00512 gGeoManager->FindNextBoundary();
00513 columnD += gGeoManager->GetStep()*node->GetMedium()->GetMaterial()->GetDensity();
00514
00516 node = gGeoManager->Step();
00517 }
00518
00521 const double *current = gGeoManager->GetCurrentPoint();
00522 length = TMath::Sqrt(TMath::Power(p2[0]-current[0], 2.)
00523 + TMath::Power(p2[1]-current[1], 2.)
00524 + TMath::Power(p2[2]-current[2], 2.));
00525 columnD += length*node->GetMedium()->GetMaterial()->GetDensity();
00526
00527 return columnD;
00528 }
00529
00530
00531
00536 void Geometry::SetDetectorSize()
00537 {
00538 double dummy;
00539
00540
00541
00542
00543
00544
00545 const TGeoVolume* vvp = gGeoManager->GetVolume("vPlaneV");
00546 const TGeoVolume* hvp = gGeoManager->GetVolume("vPlaneH");
00547
00548
00549
00550 if (vvp==0) vvp = gGeoManager->GetVolume("vPlaneVIPND");
00551 if (hvp==0) hvp = gGeoManager->GetVolume("vPlaneHIPND");
00552 if (vvp==0 && hvp==0) {
00553 vvp = gGeoManager->GetVolume("vPlaneVIPND");
00554 hvp = gGeoManager->GetVolume("vPlaneHIPND");
00555 }
00556 if (vvp==0 || hvp==0) {
00557 throw Exception(__FILE__,__LINE__,
00558 Exception::BAD_GEO_CONFIG,
00559 "Unable to find shapes to set detector size");
00560 }
00561
00562 const TGeoShape* vp = vvp->GetShape();
00563 const TGeoShape* hp = hvp->GetShape();
00564 vp->GetAxisRange(3,dummy,fDetHalfHeight);
00565 hp->GetAxisRange(3,dummy,fDetHalfWidth);
00566
00567
00568
00569 double vplanez1, vplanez2;
00570 vp->GetAxisRange(2,vplanez1,vplanez2);
00571
00572
00573
00574
00575
00576
00577
00578 double xyz1[3] = {0,0,0};
00579 double xyz2[3] = {0,0,0};
00580 this->Plane(0). Cell(0).GetCenter(xyz1);
00581 this->Plane(fPlanes.size()-1).Cell(0).GetCenter(xyz2);
00582 fDetLength = (xyz2[2]-xyz1[2])+(vplanez2-vplanez1);
00583 }
00584