Components/ATHENA/Extract_ATHENA_gdml.C

0001 #include <TFile.h>
0002 #include <TGeoManager.h>
0003 #include <TGeoNode.h>
0004 #include <TGDMLWrite.h>
0005 #include <RVersion.h>
0006
0007 #include <string>
0008
0009 int Extract_ATHENA_gdml( string subsys="all", const string outbase="athena_")
0010 {
0011   TFile * file = new TFile("detector_geometry.root","READ");
0012   gGeoManager = (TGeoManager*) file->Get("default");
0013
0014   // need to switch to root units
0015   // This flexibility was added in 6.22
0016   // And seems to have changed the default
0017 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,22,00)
0018   gGeoManager->LockDefaultUnits(false);
0019   gGeoManager->SetDefaultUnits(TGeoManager::kRootUnits);
0020 #endif
0021
0022   TGDMLWrite *writer = new TGDMLWrite;
0023   writer->SetG4Compatibility(true);
0024   writer->WriteGDMLfile ( gGeoManager, "athena_all.gdml","vg");
0025
0026   vector<string> outnames;
0027
0028   TGeoNodeMatrix* m = (TGeoNodeMatrix*) gGeoManager->GetListOfNodes()->At(0);
0029   for (int i = 0; i<m->GetNdaughters() ; ++i ){
0030     TGeoNode* n = (TGeoNode*) m->GetNodes()->At(i);
0031     string outname = outbase+n->GetName()+".gdml";
0032     if ( subsys !="all" && n->GetName()!=subsys ) continue;
0033     writer->WriteGDMLfile ( gGeoManager, n, outname.c_str(),"vg");
0034     outnames.push_back(outname);
0035   }
0036
0037   // Now need to reimport these and add a world_volume
0038   for ( auto fname : outnames ){
0039 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,22,00)
0040     gGeoManager->LockDefaultUnits(false);
0041     gGeoManager->SetDefaultUnits(TGeoManager::kRootUnits);
0042 #endif
0043     gGeoManager->Import(fname.c_str());
0044
0045     TGeoMedium *Air = gGeoManager->GetMedium("Air");
0046     if ( !Air ){
0047       TGeoElementTable *table = gGeoManager->GetElementTable();
0048       TGeoMixture *air = new TGeoMixture("air",4, 0.00120479);
0049       air->AddElement(table->GetElement(6),0.00012);
0050       air->AddElement(table->GetElement(7),0.754);
0051       air->AddElement(table->GetElement(8),0.234);
0052       air->AddElement(table->GetElement(18),0.012827);
0053       Air = new TGeoMedium("Air",0, air);
0054     }
0055
0056     TGeoVolume *top = gGeoManager->MakeBox("world_volume", Air, 3000., 3000., 10000.);
0057     top->AddNode( gGeoManager->GetTopVolume(), 1 );
0058     gGeoManager->SetTopVolume(top);
0059     gGeoManager->Export(fname.c_str(), "vg");
0060   }
0061
0062   return 0;
0063 }