Hi Stanislav,
I think the idea of being able to convert TGeometry to new TGeo is
useful (however NOT efficient - one should start directly with TGeo in
order to design an optimized geometry). Your macro does not work
correctly because your recursive method create_volumes() fills the
content of a given volume as many times as the given volume is
positioned in the geometry, resulting in huge geometries (~30 mil.
nodes) which is hard to draw even for the new system. For bigger
geometries than na49, the algorithm creates zillions of nodes and your
memory resources eventually end...
I have modified your macro so that now it seems to work. I did not
check everything in detail, but after some cleanup I think it worth
having it in /tutorials.
Thanks a lot and let me know,
Andrei.
Stanislav Nesterov wrote:
> Hi Andrei,
> I would like to convert TGeometry detector desription to TGeoManager.
> For this purpose I have written the (not yet complete) script (see
> attach.). May be it is not very efficient but I hope it does the job in
> the first approximation. I would be very grateful if you can point out
> the way to speed up it or optimize volume construction.
>
> The first problem I have encountered that resulting geometry is much
> slower then TGeometry description: I've tryed na49.root geometry file
> created by tutorial/geometry.C. And then I try to draw it.
>
> The second problem: for really big detectors CloseGeometry caused
> silent ROOT crash on 'Building caches for nodes and matrix' without any
> complain.
> ......................
> Fixing runtime shapes...
> Counting nodes...
> Voxelizing...
> Building caches for nodes and matrices...
> bash$
> .......................
>
> May be this problem is connected with not very efficient volume filling.
> So if the first issue will be solved then other will vanish.
>
> Best regards,
> Stanislav.
>
>
> ------------------------------------------------------------------------
>
> /*
> ALPHA VERSION!!!!
> DON'T USE UNLESS YOU HAVE TO!!
>
> Try converts TGeometry object int TGeoManager.
>
> To use, load libGeom.so, compile by ACLiC, load TGeometry file,
> and call
>
> convert(tgeometry,"newfilname.root")
>
>
>
> Author: S. Nesterov <Stanislav.Nesterov@cern.ch>
> */
>
> #include <TGeoManager.h>
> #include <TGeometry.h>
> #include <TNode.h>
> #include <TGeoVolume.h>
> #include <TGeoMatrix.h>
> #include <TGeoMedium.h>
> #include <TGeoMaterial.h>
> #include <TMaterial.h>
> #include <TMixture.h>
> #include <TShape.h>
> #include <TString.h>
> #include <TROOT.h>
> #include <TBRIK.h>
> #include <TTUBE.h>
> #include <TCONE.h>
> #include <TTRAP.h>
> #include <TTRD1.h>
> #include <TPGON.h>
> #include <TPCON.h>
> #include <TPARA.h>
> #include <TTUBS.h>
> #include <TCONS.h>
> #include <TELTU.h>
> #include <TGeoEltu.h>
> #include <TGeoPara.h>
> #include <TGeoPgon.h>
> #include <TGeoPcon.h>
> #include <TGeoTrd1.h>
> #include <TGeoBBox.h>
> #include <TGeoTube.h>
> #include <TGeoCone.h>
> #include <TGeoArb8.h>
> #include <Riostream.h>
> #include <TFile.h>
> TGeoManager *gGeoMan = 0;
>
> TGeoShape *makeshape(TShape*shapa){
> TString cn(shapa->ClassName());
> if(cn=="TBRIK"){
> TBRIK * bri = (TBRIK*)shapa;
> return new TGeoBBox(bri->GetDx(),bri->GetDy(),bri->GetDz());
> }
> if(cn=="TTUBE"){
> TTUBE *tb=(TTUBE*)shapa;
> return new TGeoTube(tb->GetRmin(),tb->GetRmax(),tb->GetDz());
> }
> if(cn=="TCONE"){
> TCONE *tc=(TCONE*)shapa;
> return new TGeoCone(tc->GetDz(),tc->GetRmin(),
> tc->GetRmax(),tc->GetRmin2(),tc->GetRmax2());
> }
> if(cn=="TTRAP"){
> TTRAP *trp=(TTRAP*)shapa;
> return new TGeoTrap(trp->GetDz(),trp->GetDx(),trp->GetDy(),
> trp->GetH1(),trp->GetBl1(),trp->GetTl1(),
> trp->GetAlpha1(),trp->GetH2(),trp->GetBl2(),
> trp->GetTl2(),trp->GetAlpha2());
>
> }
> if(cn=="TTRD1"){
> TTRD1 *trd=(TTRD1*)shapa;
> return new TGeoTrd1(trd->GetDx(),trd->GetDx2(),trd->GetDy(),trd->GetDz());
> }
> if(cn=="TPGON"){
> TPGON *pgn=(TPGON*)shapa;
> Int_t nz = pgn->GetNz();
> TGeoPgon*gpg= new TGeoPgon(pgn->GetPhi1(),pgn->GetDhi1(),
> pgn->GetNdiv(),nz);
> Float_t *rmin = pgn->GetRmin();
> Float_t *rmax = pgn->GetRmax();
> Float_t *dz = pgn->GetDz();
> for(int i=0;i<nz;i++){
> gpg->DefineSection(i,dz[i],rmin[i],rmax[i]);
> }
> return gpg;
> }
> if(cn=="TPCON"){
> TPCON *pgn=(TPCON*)shapa;
> Int_t nz = pgn->GetNz();
> TGeoPcon *gpg = new TGeoPcon(pgn->GetPhi1(),pgn->GetDhi1(),nz);
> Float_t *rmin = pgn->GetRmin();
> Float_t *rmax = pgn->GetRmax();
> Float_t *dz = pgn->GetDz();
> for(int i=0;i<nz;i++){
> gpg->DefineSection(i,dz[i],rmin[i],rmax[i]);
> }
> return gpg;
> }
> if(cn=="TPARA"){
> TPARA *par=(TPARA*)shapa;
> return new TGeoPara(par->GetDx(),par->GetDy(),par->GetDz(),
> par->GetAlpha(),par->GetTheta(),par->GetPhi());
> }
> if(cn=="TTUBS"){
> TTUBS *tbs=(TTUBS*)shapa;
> return new TGeoTubeSeg(tbs->GetRmin(),tbs->GetRmax(),tbs->GetDz(),
> tbs->GetPhi1(),tbs->GetPhi2());
> }
> if(cn=="TCONS"){
> TCONS *tcs=(TCONS*)shapa;
> return new TGeoConeSeg(tcs->GetDz(),tcs->GetRmin(),tcs->GetRmax(),
> tcs->GetRmin2(),tcs->GetRmax2(),
> tcs->GetPhi1(),tcs->GetPhi2());
> }
>
> if(cn=="TELTU"){
> TELTU *tut=(TELTU*)shapa;
> return new TGeoEltu(tut->GetRmax(),tut->GetRmax()*tut->GetAspectRatio(),
> tut->GetDz());
> }/*
> if(cn=="THYPE"){
> THYPE *hyp=(THYPE*)shapa;
> }
>
> if(cn=="TGTRA"){
> TGTRA *gtr=(TGTRA*)shapa;
> }
> if(cn=="TTRD2"){
> TTRD2 *td2=(TTRD2*)shapa;
> }
>
> if(cn=="TXTRU"){
> TXTRU *xtr=(TXTRU*)shapa;
> }
> */ cerr<<"Non-implemented shape: "<<cn<<endl;
> return 0;
> if(cn=="TCTUB"){
> // TCTUB *tct=(TCTUB*)shapa;
>
> }
> }
> void create_volumes(TNode *nd=0,TGeoVolume* volm=0){
> if(!nd ) return;
> if (!volm) return;
> TIter iter(nd->GetListOfNodes());
>
> cout << ":"<<nd->GetName()<<"(";
>
> TNode *next = 0;
> while((next = (TNode*)iter())) {
>
> // Search for translation and rotation transformations
> TRotMatrix *mtr = next->GetMatrix();
> TGeoRotation *grot = 0;
> TGeoMatrix *gmat = 0;
> TGeoTranslation *gtran = new TGeoTranslation(next->GetX(),next->GetY(),
> next->GetZ());
> if(mtr->GetType()){ // if the matrix is not identity
> grot = new TGeoRotation(mtr->GetName());
> grot->SetMatrix(mtr->GetMatrix());
>
>
> gmat = new TGeoCombiTrans(next->GetX(),next->GetY(),next->GetZ(),
> grot);
> // gmat = new TGeoCombiTrans(*gtran,*grot);
> } else
> gmat = gtran;
> TGeoVolume *curvol = gGeoMan->GetVolume(next->GetShape()->GetName());
> if(!curvol){
> cerr<<"No such shape!!: "<<next->GetShape()->GetName()<<endl;
> }
> if(next->GetVisibility()>0) curvol->SetVisibility(kTRUE);
> else curvol->SetVisibility(kFALSE);
> create_volumes(next,curvol);
> curvol->SetNumber(curvol->GetNumber()+1);
> volm->AddNode(curvol,curvol->GetNumber(),gmat);
> }
> cout<<")";
> }
>
> void convert(TGeometry *geom=0,char *name="tmp.root"){
> // gSystem->Load("libGeom.so");
> TFile stmp(name,"recreate");
> gGeoMan = new TGeoManager("geoman","Geomanager");
> if(!geom) return;
> cout<<"Defines media"<<endl;
> TIter maxt(geom->GetListOfMaterials());
> TMaterial *mata = 0;
> while ((mata = (TMaterial*)maxt())){
> TGeoMaterial *geomat = 0;
> if(mata->ClassName() != "TGeoMixture") {
> geomat = new TGeoMaterial(mata->GetName(),mata->GetA(),
> mata->GetZ(),mata->GetDensity(),
> mata->GetRadLength(),
> mata->GetInterLength());
> } else {
> TMixture *matax = (TMixture*)mata;
> geomat = new TGeoMixture(matax->GetName(),matax->GetNmixt());
>
> for(int i=0;i<matax->GetNmixt();i++)
>
> ((TGeoMixture*)geomat)->DefineElement(i,matax->GetAmixt()[i],matax->GetZmixt()[i],
> matax->GetWmixt()[i]);
> }
> TGeoMedium *geomed = new TGeoMedium(mata->GetName(),mata->GetNumber(),
> geomat);
> geomed = 0;
>
> }
> cout<<"Defines shapes"<<endl;
> TIter voxt(geom->GetListOfShapes());
> TShape *shp = 0;
> while ((shp = (TShape*)voxt())){
> TGeoMedium *gemed = gGeoMan->GetMedium(shp->GetMaterial()->GetName());
> TGeoVolume *ngevol = new TGeoVolume(shp->GetName(),makeshape(shp),
> gemed);
> ngevol->SetLineColor(shp->GetLineColor());
> }
> cout<<"Defines volumes"<<endl;
> TIter noxt(geom->GetListOfNodes());
> TNode *nd = (TNode*)noxt();
> TGeoMaterial *va = new TGeoMaterial("Vacuum",1,0,0,0);
> TGeoMedium *vac = new TGeoMedium("Vacuum",1,va);
> TGeoShape *shap = makeshape(nd->GetShape());
> if(!shap) {
> cout<<"Shape is not implemented"<<endl;
> return;
> }
> cout<<"Create hierarchy"<<endl;
> TGeoVolume *top = new TGeoVolume("TOP",shap,vac);
> top->SetVisibility(kFALSE);
> gGeoMan->SetTopVolume(top);
> /* // just debuggin lines
> TGeoVolume*gl = gGeoMan->GetVolume("SLGC");
> top->AddNode(gl,1, gGeoIdentity); */
>
> // Main recursive procedure to fill geometry volume
> // VERY CPU consuming procedure
> create_volumes(nd,top);
>
> // Closing geometry - quite lengthy
> gGeoMan->CloseGeometry();
>
> gGeoMan->SetVisLevel(4);
>
> gGeoMan->Write();
> stmp.Close();
> }
>
/*
ALPHA VERSION!!!!
DON'T USE UNLESS YOU HAVE TO!!
Try converts TGeometry object int TGeoManager.
To use, load libGeom.so, compile by ACLiC, load TGeometry file,
and call
convert(tgeometry,"newfilname.root")
Author: S. Nesterov <Stanislav.Nesterov@cern.ch>
*/
#include <TGeoManager.h>
#include <TGeometry.h>
#include <TNode.h>
#include <TGeoVolume.h>
#include <TGeoMatrix.h>
#include <TGeoMedium.h>
#include <TGeoMaterial.h>
#include <TMaterial.h>
#include <TMixture.h>
#include <TShape.h>
#include <TString.h>
#include <TROOT.h>
#include <TBRIK.h>
#include <TTUBE.h>
#include <TCONE.h>
#include <TTRAP.h>
#include <TTRD1.h>
#include <TPGON.h>
#include <TPCON.h>
#include <TPARA.h>
#include <TTUBS.h>
#include <TCONS.h>
#include <TELTU.h>
#include <TGeoEltu.h>
#include <TGeoPara.h>
#include <TGeoPgon.h>
#include <TGeoPcon.h>
#include <TGeoTrd1.h>
#include <TGeoBBox.h>
#include <TGeoTube.h>
#include <TGeoCone.h>
#include <TGeoArb8.h>
#include <Riostream.h>
#include <TFile.h>
TGeoManager *gGeoMan = 0;
TGeoShape *makeshape(TShape*shapa){
TString cn(shapa->ClassName());
if(cn=="TBRIK"){
TBRIK * bri = (TBRIK*)shapa;
return new TGeoBBox(bri->GetDx(),bri->GetDy(),bri->GetDz());
}
if(cn=="TTUBE"){
TTUBE *tb=(TTUBE*)shapa;
return new TGeoTube(tb->GetRmin(),tb->GetRmax(),tb->GetDz());
}
if(cn=="TCONE"){
TCONE *tc=(TCONE*)shapa;
return new TGeoCone(tc->GetDz(),tc->GetRmin(),
tc->GetRmax(),tc->GetRmin2(),tc->GetRmax2());
}
if(cn=="TTRAP"){
TTRAP *trp=(TTRAP*)shapa;
return new TGeoTrap(trp->GetDz(),trp->GetDx(),trp->GetDy(),
trp->GetH1(),trp->GetBl1(),trp->GetTl1(),
trp->GetAlpha1(),trp->GetH2(),trp->GetBl2(),
trp->GetTl2(),trp->GetAlpha2());
}
if(cn=="TTRD1"){
TTRD1 *trd=(TTRD1*)shapa;
return new TGeoTrd1(trd->GetDx(),trd->GetDx2(),trd->GetDy(),trd->GetDz());
}
if(cn=="TPGON"){
TPGON *pgn=(TPGON*)shapa;
Int_t nz = pgn->GetNz();
TGeoPgon*gpg= new TGeoPgon(pgn->GetPhi1(),pgn->GetDhi1(),
pgn->GetNdiv(),nz);
Float_t *rmin = pgn->GetRmin();
Float_t *rmax = pgn->GetRmax();
Float_t *dz = pgn->GetDz();
for(int i=0;i<nz;i++){
gpg->DefineSection(i,dz[i],rmin[i],rmax[i]);
}
return gpg;
}
if(cn=="TPCON"){
TPCON *pgn=(TPCON*)shapa;
Int_t nz = pgn->GetNz();
TGeoPcon *gpg = new TGeoPcon(pgn->GetPhi1(),pgn->GetDhi1(),nz);
Float_t *rmin = pgn->GetRmin();
Float_t *rmax = pgn->GetRmax();
Float_t *dz = pgn->GetDz();
for(int i=0;i<nz;i++){
gpg->DefineSection(i,dz[i],rmin[i],rmax[i]);
}
return gpg;
}
if(cn=="TPARA"){
TPARA *par=(TPARA*)shapa;
return new TGeoPara(par->GetDx(),par->GetDy(),par->GetDz(),
par->GetAlpha(),par->GetTheta(),par->GetPhi());
}
if(cn=="TTUBS"){
TTUBS *tbs=(TTUBS*)shapa;
return new TGeoTubeSeg(tbs->GetRmin(),tbs->GetRmax(),tbs->GetDz(),
tbs->GetPhi1(),tbs->GetPhi2());
}
if(cn=="TCONS"){
TCONS *tcs=(TCONS*)shapa;
return new TGeoConeSeg(tcs->GetDz(),tcs->GetRmin(),tcs->GetRmax(),
tcs->GetRmin2(),tcs->GetRmax2(),
tcs->GetPhi1(),tcs->GetPhi2());
}
if(cn=="TELTU"){
TELTU *tut=(TELTU*)shapa;
return new TGeoEltu(tut->GetRmax(),tut->GetRmax()*tut->GetAspectRatio(),
tut->GetDz());
}/*
if(cn=="THYPE"){
THYPE *hyp=(THYPE*)shapa;
}
if(cn=="TGTRA"){
TGTRA *gtr=(TGTRA*)shapa;
}
if(cn=="TTRD2"){
TTRD2 *td2=(TTRD2*)shapa;
}
if(cn=="TXTRU"){
TXTRU *xtr=(TXTRU*)shapa;
}
*/ cerr<<"Non-implemented shape: "<<cn<<endl;
return 0;
if(cn=="TCTUB"){
// TCTUB *tct=(TCTUB*)shapa;
}
}
void create_volumes1(TNode *node=0, TGeoVolume *volm=0){
if (!node) return;
if (!volm) volm = gGeoManager->GetVolume(node->GetShape()->GetName());
if (!volm) {
cerr << "No volume " << node->GetShape()->GetName()<<endl;
return;
}
if (volm->GetNdaughters()) return;
TList *dlist = node->GetListOfNodes();
if (!dlist) return;
Int_t nd = dlist->GetSize();
if (!nd) return;
Int_t id;
TNode *daughter;
TGeoVolume *vold;
for (id=0; id<nd; id++) {
daughter = (TNode*)dlist->At(id);
if (!daughter) return; // normally useless
vold = gGeoManager->GetVolume(daughter->GetShape()->GetName());
// printf("VOL:%s (%i) adding node: %s_%i\n", volm->GetName(), volm->GetNumber(),vold->GetName(), id);
if (!vold) {
cerr << "No volume " << daughter->GetShape()->GetName()<<endl;
return;
}
TRotMatrix *mtr = daughter->GetMatrix();
TGeoRotation *grot = 0;
TGeoMatrix *gmat = 0;
Double_t x,y,z;
x = daughter->GetX();
y = daughter->GetY();
z = daughter->GetZ();
Bool_t is_translated = ((x==0) && (y==0) && (z==0))?kFALSE:kTRUE;
Bool_t is_rotated = (mtr->GetType())?kTRUE:kFALSE;
// if ((x==0) && (y==0) && (z==0)) gtran = new TGeoTranslation(x,y,z);
if(is_rotated){ // if the matrix is not identity
grot = new TGeoRotation(mtr->GetName());
grot->SetMatrix(mtr->GetMatrix());
if (is_translated)
gmat = new TGeoCombiTrans(x,y,z,grot);
else
gmat = grot;
} else {
if (is_translated)
gmat = new TGeoTranslation(x,y,z);
else
gmat = gGeoIdentity;
}
volm->AddNode(vold, id, gmat);
create_volumes1(daughter, vold);
}
}
void create_volumes(TNode *nd=0,TGeoVolume* volm=0){
if(!nd ) return;
if (!volm) return;
TIter iter(nd->GetListOfNodes());
cout << ":"<<nd->GetName()<<"(";
TNode *next = 0;
while((next = (TNode*)iter())) {
// Search for translation and rotation transformations
TRotMatrix *mtr = next->GetMatrix();
TGeoRotation *grot = 0;
TGeoMatrix *gmat = 0;
TGeoTranslation *gtran = new TGeoTranslation(next->GetX(),next->GetY(),
next->GetZ());
if(mtr->GetType()){ // if the matrix is not identity
grot = new TGeoRotation(mtr->GetName());
grot->SetMatrix(mtr->GetMatrix());
gmat = new TGeoCombiTrans(next->GetX(),next->GetY(),next->GetZ(),
grot);
// gmat = new TGeoCombiTrans(*gtran,*grot);
} else
gmat = gtran;
TGeoVolume *curvol = gGeoMan->GetVolume(next->GetShape()->GetName());
if(!curvol){
cerr<<"No such shape!!: "<<next->GetShape()->GetName()<<endl;
}
if(next->GetVisibility()>0) curvol->SetVisibility(kTRUE);
else curvol->SetVisibility(kFALSE);
create_volumes(next,curvol);
curvol->SetNumber(curvol->GetNumber()+1);
volm->AddNode(curvol,curvol->GetNumber(),gmat);
}
cout<<")";
}
void convert(TGeometry *geom=0,char *name="tmp.root"){
// gSystem->Load("libGeom.so");
// TFile stmp(name,"recreate");
gGeoMan = new TGeoManager("geoman","Geomanager");
if(!geom) return;
cout<<"Defines media"<<endl;
TIter maxt(geom->GetListOfMaterials());
TMaterial *mata = 0;
while ((mata = (TMaterial*)maxt())){
TGeoMaterial *geomat = 0;
if(mata->ClassName() != "TGeoMixture") {
geomat = new TGeoMaterial(mata->GetName(),mata->GetA(),
mata->GetZ(),mata->GetDensity(),
mata->GetRadLength(),
mata->GetInterLength());
} else {
TMixture *matax = (TMixture*)mata;
geomat = new TGeoMixture(matax->GetName(),matax->GetNmixt());
for(int i=0;i<matax->GetNmixt();i++)
((TGeoMixture*)geomat)->DefineElement(i,matax->GetAmixt()[i],matax->GetZmixt()[i],
matax->GetWmixt()[i]);
}
TGeoMedium *geomed = new TGeoMedium(mata->GetName(),mata->GetNumber(),
geomat);
geomed = 0;
}
cout<<"Defines shapes"<<endl;
TIter voxt(geom->GetListOfShapes());
TShape *shp = 0;
while ((shp = (TShape*)voxt())){
TGeoMedium *gemed = gGeoMan->GetMedium(shp->GetMaterial()->GetName());
TGeoVolume *ngevol = new TGeoVolume(shp->GetName(),makeshape(shp),
gemed);
ngevol->SetLineColor(shp->GetLineColor());
}
cout<<"Defines volumes"<<endl;
TIter noxt(geom->GetListOfNodes());
TNode *nd = (TNode*)noxt();
cout<<"Create hierarchy"<<endl;
TGeoVolume *top = gGeoManager->GetVolume(nd->GetShape()->GetName());
top->SetVisibility(kFALSE);
gGeoMan->SetTopVolume(top);
/* // just debuggin lines
TGeoVolume*gl = gGeoMan->GetVolume("SLGC");
top->AddNode(gl,1, gGeoIdentity); */
// Main recursive procedure to fill geometry volume
// VERY CPU consuming procedure
// --- yes, if a volume was placed in several others, you were adding the same
// content several times, resulting in huge geometries (~30 mil nodes!!!)
// create_volumes(nd,top);
create_volumes1(nd,top);
// Closing geometry - quite lengthy
gGeoMan->CloseGeometry();
gGeoMan->SetVisLevel(4);
// TGeoManager::Write() does not work for the time being (see Export/Import)
gGeoMan->Export(name);
// stmp.Close();
}
This archive was generated by hypermail 2b29 : Thu Jan 01 2004 - 17:50:13 MET