Skip to content
Snippets Groups Projects
Commit 1881ba46 authored by Markus Frank's avatar Markus Frank
Browse files

DD4hep Geant4 geometry converter

parent 2f214b5b
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
DDG4/CMakeLists.txt 0 → 100644
+ 18
− 0
View file @ 1881ba46
#---Find Geant4-------------------------------------------------------------------
find_package(Geant4 REQUIRED)
if(NOT Geant4_clhep_FOUND)
find_package(CLHEP REQUIRED)
set(Geant4_INCLUDE_DIRS ${Geant4_INCLUDE_DIRS} ${CLHEP_INCLUDE_DIRS})
set(Geant4_LIBRARIES ${Geant4_LIBRARIES} ${CLHEP_LIBRARIES})
endif()
#---Includedirs-------------------------------------------------------------------
include_directories(${CMAKE_SOURCE_DIR}/DDCore/include
${CMAKE_CURRENT_SOURCE_DIR}/include
${ROOT_INCLUDE_DIR}
${Geant4_INCLUDE_DIRS})
#---Add Library-------------------------------------------------------------------
file(GLOB sources src/*.cpp)
add_library(DD4hepG4 SHARED ${sources})
target_link_libraries(DD4hepG4 DD4hepCore ${ROOT_LIBRARIES} ${Geant4_LIBRARIES})
DDG4/include/Geant4Converter.h 0 → 100644
+ 125
− 0
View file @ 1881ba46
// $Id:$
//====================================================================
// AIDA Detector description implementation
//--------------------------------------------------------------------
//
// Author : M.Frank
//
//====================================================================
#ifndef DD4HEP_GEANT4CONVERTER_H
#define DD4HEP_GEANT4CONVERTER_H
#include "DD4hep/GeoHandler.h"
#include "DD4hep/LCDD.h"
#include <set>
#include <map>
#include <vector>
class TGeoVolume;
class TGeoNode;
#if 0
struct G4VAny { virtual ~G4VAny() {} };
template <class T> class G4AnyThing : public T {
public:
G4AnyThing() : T() {}
template <class A> G4AnyThing(A) : T() {}
template <class A,class B> G4AnyThing(A,B) : T() {}
template <class A,class B,class C> G4AnyThing(A,B,C) : T() {}
template <class A, class B, class C, class D> G4AnyThing(A,B,C,D) : T() {}
template <class A, class B, class C, class D, class E> G4AnyThing(A,B,C,D,E) : T() {}
template <class A, class B, class C, class D, class E, class F> G4AnyThing(A,B,C,D,E,F) : T() {}
template <class A, class B, class C, class D, class E, class F, class G> G4AnyThing(A,B,C,D,E,F,G) : T() {}
template <class A, class B, class C, class D, class E, class F, class G, class H> G4AnyThing(A,B,C,D,E,F,G,H) : T() {}
template <class A, class B, class C, class D, class E, class F, class G, class H, class I> G4AnyThing(A,B,C,D,E,F,G,H,I) : T() {}
void AddIsotope(void*, double) {}
void AddElement(void*, double) {}
static G4AnyThing<G4VAny>* GetIsotope(const std::string&, bool) { return 0; }
static G4AnyThing<G4VAny>* GetElement(const std::string&, bool) { return 0; }
static G4AnyThing<G4VAny>* GetMaterial(const std::string&, bool) { return 0; }
int GetNelements() const { return 0; }
double GetDensity() const { return 0.0; }
};
typedef G4VAny G4VSolid;
typedef G4AnyThing<G4VSolid> G4Box;
typedef G4AnyThing<G4VSolid> G4Tubs;
typedef G4AnyThing<G4VSolid> G4Trd;
typedef G4AnyThing<G4VSolid> G4Paraboloid;
typedef G4AnyThing<G4VSolid> G4Polycone;
typedef G4AnyThing<G4VSolid> G4Sphere;
typedef G4AnyThing<G4VSolid> G4Torus;
typedef G4AnyThing<G4VSolid> G4UnionSolid;
typedef G4AnyThing<G4VSolid> G4SubtractionSolid;
typedef G4AnyThing<G4VSolid> G4IntersectionSolid;
typedef G4AnyThing<G4VAny> G4LogicalVolume;
typedef G4AnyThing<G4VAny> G4Material;
typedef G4AnyThing<G4VAny> G4Element;
typedef G4AnyThing<G4VAny> G4Isotope;
typedef G4AnyThing<G4VAny> G4Transform3D;
typedef G4AnyThing<G4VAny> G4ThreeVector;
typedef G4AnyThing<G4VAny> G4PVPlacement;
#else
#include "G4Element.hh"
#include "G4Box.hh"
#include "G4Tubs.hh"
#include "G4Trd.hh"
#include "G4Paraboloid.hh"
#include "G4Polycone.hh"
#include "G4Sphere.hh"
#include "G4Torus.hh"
#include "G4UnionSolid.hh"
#include "G4SubtractionSolid.hh"
#include "G4IntersectionSolid.hh"
#include "G4LogicalVolume.hh"
#include "G4Material.hh"
#include "G4Element.hh"
#include "G4Isotope.hh"
#include "G4Transform3D.hh"
#include "G4ThreeVector.hh"
#include "G4PVPlacement.hh"
#endif
namespace DD4hep {
namespace Geometry {
struct Geant4Converter : public GeoHandler {
struct G4GeometryInfo : public GeometryInfo {
std::map<const TGeoElement*,G4Element*> g4Elements;
std::map<const TGeoMedium*, G4Material*> g4Materials;
std::map<const TGeoShape*, G4VSolid*> g4Solids;
std::map<const TGeoVolume*, G4LogicalVolume*> g4Volumes;
std::map<const TGeoNode*, G4PVPlacement*> g4Placements;
};
G4GeometryInfo* m_dataPtr;
G4GeometryInfo& data() const { return *m_dataPtr; }
/// Constructor
Geant4Converter() {}
/// Standard destructor
virtual ~Geant4Converter() {}
/// Create geometry dump
void create(DetElement top);
/// Dump material in GDML format to output stream
virtual void* handleMaterial(const std::string& name, const TGeoMedium* medium) const;
/// Dump element in GDML format to output stream
virtual void* handleElement(const std::string& name, const TGeoElement* element) const;
/// Dump solid in GDML format to output stream
virtual void* handleSolid(const std::string& name, const TGeoShape* volume) const;
/// Dump logical volume in GDML format to output stream
virtual void* handleVolume(const std::string& name, const TGeoVolume* volume) const;
/// Dump volume placement in GDML format to output stream
virtual void* handlePlacement(const std::string& name, const TGeoNode* node) const;
};
} // End namespace Geometry
} // End namespace DD4hep
#endif // DD4HEP_GEANT4CONVERTER_H
DDG4/src/Geant4Converter.cpp 0 → 100644
+ 313
− 0
View file @ 1881ba46
// $Id:$
//====================================================================
// AIDA Detector description implementation for LCD
//--------------------------------------------------------------------
//
// Author : M.Frank
//
//====================================================================
#include "DD4hep/LCDD.h"
#include "Geant4Converter.h"
// ROOT includes
#include "TROOT.h"
#include "TColor.h"
#include "TGeoShape.h"
#include "TGeoCone.h"
#include "TGeoParaboloid.h"
#include "TGeoPcon.h"
#include "TGeoPgon.h"
#include "TGeoSphere.h"
#include "TGeoTorus.h"
#include "TGeoTube.h"
#include "TGeoTrd1.h"
#include "TGeoTrd2.h"
#include "TGeoArb8.h"
#include "TGeoMatrix.h"
#include "TGeoBoolNode.h"
#include "TGeoCompositeShape.h"
#include "TGeoNode.h"
#include "TClass.h"
#include "TMath.h"
#include <iostream>
using namespace DD4hep::Geometry;
using namespace DD4hep;
using namespace std;
namespace {
static string indent = "";
struct MyTransform3D : public G4Transform3D {
MyTransform3D(double XX, double XY, double XZ, double DX,
double YX, double YY, double YZ, double DY,
double ZX, double ZY, double ZZ, double DZ) : G4Transform3D()
{}
// { G4Transform3D::setTransform(XX,XY,XZ,DX,YX,YY,YZ,DY,ZX,ZY,ZZ,DZ); }
};
}
/// Dump element in GDML format to output stream
void* Geant4Converter::handleElement(const string& name, const TGeoElement* element) const {
G4Element* g4e = data().g4Elements[element];
if ( !g4e ) {
g4e = G4Element::GetElement(name,false);
if ( !g4e ) {
if ( element->GetNisotopes() > 1 ) {
g4e = new G4Element(name,element->GetTitle(),element->GetNisotopes());
for(int i=0, n=element->GetNisotopes(); i<n; ++i) {
TGeoIsotope* iso = element->GetIsotope(i);
G4Isotope* g4iso = G4Isotope::GetIsotope(iso->GetName(),false);
if ( !g4iso ) {
g4iso = new G4Isotope(iso->GetName(),iso->GetZ(),iso->GetN(),iso->GetA());
}
g4e->AddIsotope(g4iso,element->GetRelativeAbundance(i));
}
}
else {
g4e = new G4Element(name,element->GetTitle(),element->A(),element->Z());
}
}
data().g4Elements[element] = g4e;
}
return g4e;
}
/// Dump material in GDML format to output stream
void* Geant4Converter::handleMaterial(const std::string& name, const TGeoMedium* medium) const {
G4Material* mat = data().g4Materials[medium];
if ( !mat ) {
mat = G4Material::GetMaterial(name,false);
if ( !mat ) {
TGeoMaterial* m = medium->GetMaterial();
if ( m->IsMixture() ) {
TGeoMixture* mix = (TGeoMixture*)m;
mat = new G4Material(name,mix->GetDensity(),mix->GetNelements());
for(int i=0, n=mix->GetNelements(); i<n; ++i) {
TGeoElement* e = mix->GetElement(i);
G4Element* g4e = (G4Element*)handleElement(e->GetName(),e);
mat->AddElement(g4e,(mix->GetWmixt())[i]);
}
}
else {
mat = new G4Material(name,m->GetZ(),m->GetA(),m->GetDensity());
}
}
data().g4Materials[medium] = mat;
}
return mat;
}
/// Dump solid in GDML format to output stream
void* Geant4Converter::handleSolid(const string& name, const TGeoShape* shape) const {
G4VSolid* solid = data().g4Solids[shape];
if ( !solid && shape ) {
if ( shape->IsA() == TGeoBBox::Class() ) {
const TGeoBBox* s = (const TGeoBBox*)shape;
solid = new G4Box(name,s->GetDX()*CM_2_MM,s->GetDY()*CM_2_MM,s->GetDZ()*CM_2_MM);
}
else if ( shape->IsA() == TGeoTube::Class() ) {
const TGeoTube* s = (const TGeoTube*)shape;
solid = new G4Tubs(name,s->GetRmin()*CM_2_MM,s->GetRmax()*CM_2_MM,s->GetDz()*CM_2_MM,0,2.*M_PI);
}
else if ( shape->IsA() == TGeoTubeSeg::Class() ) {
const TGeoTubeSeg* s = (const TGeoTubeSeg*)shape;
solid = new G4Tubs(name,s->GetRmin()*CM_2_MM,s->GetRmax()*CM_2_MM,s->GetDz()*CM_2_MM,s->GetPhi1()*DEGREE_2_RAD,s->GetPhi2()*DEGREE_2_RAD);
}
else if ( shape->IsA() == TGeoTrd1::Class() ) {
const TGeoTrd1* s = (const TGeoTrd1*)shape;
solid = new G4Trd(name,s->GetDx1()*CM_2_MM,s->GetDx2()*CM_2_MM,s->GetDy()*CM_2_MM,s->GetDy()*CM_2_MM,s->GetDz()*CM_2_MM);
}
else if ( shape->IsA() == TGeoTrd2::Class() ) {
const TGeoTrd2* s = (const TGeoTrd2*)shape;
solid = new G4Trd(name,s->GetDx1()*CM_2_MM,s->GetDx2()*CM_2_MM,s->GetDy1()*CM_2_MM,s->GetDy2()*CM_2_MM,s->GetDz()*CM_2_MM);
}
else if ( shape->IsA() == TGeoPgon::Class() ) {
#if 0
const TGeoPgon* s = (const TGeoPgon*)shape;
#endif
const TGeoPcon* s = (const TGeoPcon*)shape;
vector<double> rmin, rmax, z;
for( size_t i=0; i<s->GetNz(); ++i ) {
rmin.push_back(s->GetRmin(i)*CM_2_MM);
rmax.push_back(s->GetRmax(i)*CM_2_MM);
z.push_back(s->GetZ(i)*CM_2_MM);
}
solid = new G4Polycone(name,s->GetPhi1()*DEGREE_2_RAD,(s->GetDphi()+s->GetPhi1())*DEGREE_2_RAD,s->GetNz(),
&rmin[0], &rmax[0], &z[0]);
}
else if ( shape->IsA() == TGeoPcon::Class() ) {
const TGeoPcon* s = (const TGeoPcon*)shape;
vector<double> rmin, rmax, z;
for( size_t i=0; i<s->GetNz(); ++i ) {
rmin.push_back(s->GetRmin(i)*CM_2_MM);
rmax.push_back(s->GetRmax(i)*CM_2_MM);
z.push_back(s->GetZ(i)*CM_2_MM);
}
solid = new G4Polycone(name,s->GetPhi1()*DEGREE_2_RAD,(s->GetDphi()+s->GetPhi1())*DEGREE_2_RAD,s->GetNz(),
&rmin[0], &rmax[0], &z[0]);
}
else if ( shape->IsA() == TGeoParaboloid::Class() ) {
const TGeoParaboloid* s = (const TGeoParaboloid*)shape;
solid = new G4Paraboloid(name,s->GetDz()*CM_2_MM,s->GetRlo()*CM_2_MM,s->GetRhi()*CM_2_MM);
}
else if ( shape->IsA() == TGeoSphere::Class() ) {
const TGeoSphere* s = (const TGeoSphere*)shape;
solid = new G4Sphere(name,s->GetRmin()*CM_2_MM,s->GetRmax()*CM_2_MM,
s->GetPhi1()*DEGREE_2_RAD,s->GetPhi2()*DEGREE_2_RAD,
s->GetTheta1()*DEGREE_2_RAD,s->GetTheta2()*DEGREE_2_RAD);
}
else if ( shape->IsA() == TGeoTorus::Class() ) {
const TGeoTorus* s = (const TGeoTorus*)shape;
solid = new G4Torus(name,s->GetRmin()*CM_2_MM,s->GetRmax()*CM_2_MM, s->GetR()*CM_2_MM,
s->GetPhi1()*DEGREE_2_RAD,s->GetDphi()*DEGREE_2_RAD);
}
else if ( shape->IsA() == TGeoCompositeShape::Class() ) {
const TGeoCompositeShape* s = (const TGeoCompositeShape*)shape;
const TGeoBoolNode* boolean = s->GetBoolNode();
TGeoBoolNode::EGeoBoolType oper = boolean->GetBooleanOperator();
TGeoMatrix* m = boolean->GetRightMatrix();
G4VSolid* left = (G4VSolid*)handleSolid(name+"_left", boolean->GetLeftShape());
G4VSolid* right = (G4VSolid*)handleSolid(name+"_right",boolean->GetRightShape());
const Double_t *t = m->GetTranslation();
const Double_t *r = m->GetRotationMatrix();
if ( !left ) {
throw runtime_error("G4Converter: No left Geant4 Solid present for composite shape:"+name);
}
if ( !right ) {
throw runtime_error("G4Converter: No right Geant4 Solid present for composite shape:"+name);
}
if ( m->IsRotation() ) {
MyTransform3D transform(r[0],r[1],r[2],t[0],
r[3],r[4],r[5],t[1],
r[6],r[7],r[8],t[3]);
if ( oper == TGeoBoolNode::kGeoSubtraction )
solid = new G4SubtractionSolid(name,left,right,transform);
else if ( oper == TGeoBoolNode::kGeoUnion )
solid = new G4UnionSolid(name,left,right,transform);
else if ( oper == TGeoBoolNode::kGeoIntersection )
solid = new G4IntersectionSolid(name,left,right,transform);
}
else {
G4ThreeVector transform(t[0],t[1],t[2]);
if ( oper == TGeoBoolNode::kGeoSubtraction )
solid = new G4SubtractionSolid(name,left,right,0,transform);
else if ( oper == TGeoBoolNode::kGeoUnion )
solid = new G4UnionSolid(name,left,right,0,transform);
else if ( oper == TGeoBoolNode::kGeoIntersection )
solid = new G4IntersectionSolid(name,left,right,0,transform);
}
}
if ( !solid ) {
string err = "Failed to handle unknown solid shape:" +
name + " of type " + string(shape->IsA()->GetName());
throw runtime_error(err);
}
data().g4Solids[shape] = solid;
}
return solid;
}
/// Dump logical volume in GDML format to output stream
void* Geant4Converter::handleVolume(const string& name, const TGeoVolume* volume) const {
const TGeoVolume* v = volume;
G4LogicalVolume* vol = data().g4Volumes[v];
if ( !vol ) {
string n = v->GetName();
TGeoMedium* m = v->GetMedium();
TGeoShape* s = v->GetShape();
G4VSolid* solid = (G4VSolid*)handleSolid(s->GetName(),s);
G4Material* medium = (G4Material*)handleMaterial(m->GetName(),m);
if ( !solid ) {
throw runtime_error("G4Converter: No Geant4 Solid present for volume:"+n);
}
if ( !medium ) {
throw runtime_error("G4Converter: No Geant4 material present for volume:"+n);
}
vol = new G4LogicalVolume(solid, medium, n);
data().g4Volumes[v] = vol;
}
return vol;
}
/// Dump volume placement in GDML format to output stream
void* Geant4Converter::handlePlacement(const std::string& name, const TGeoNode* node) const {
G4PVPlacement* g4pv = data().g4Placements[node];
if ( !g4pv ) {
TGeoMatrix* trafo = node->GetMatrix();
if ( trafo ) {
const Double_t* trans = trafo->GetTranslation();
const Double_t* rot = trafo->GetRotationMatrix();
int copy = node->GetNumber();
G4LogicalVolume* g4vol = data().g4Volumes[node->GetVolume()];
G4LogicalVolume* g4mot = data().g4Volumes[node->GetMotherVolume()];
if ( trafo->IsRotation() ) {
MyTransform3D transform(rot[0],rot[1],rot[2],trans[0],
rot[3],rot[4],rot[5],trans[1],
rot[6],rot[7],rot[8],trans[3]);
g4pv = new G4PVPlacement(transform, // no rotation
g4vol, // its logical volume
name, // its name
g4mot, // its mother (logical) volume
false, // no boolean operations
copy); // its copy number
}
else {
G4ThreeVector pos(trans[0],trans[1],trans[2]);
g4pv = new G4PVPlacement(0, // no rotation
pos, // translation position
g4vol, // its logical volume
name, // its name
g4mot, // its mother (logical) volume
false, // no boolean operations
copy); // its copy number
}
data().g4Placements[node] = g4pv;
// cout << "Created volume placement:" << name << " No:" << copy << endl;
}
}
else {
cout << "Attempt to DOUBLE-place physical volume:" << name << " No:" << node->GetNumber() << endl;
}
return g4pv;
}
template <typename O, typename C, typename F> void handle(const O* o, const C& c, F pmf) {
for(typename C::const_iterator i=c.begin(); i != c.end(); ++i) {
(o->*pmf)((*i).first, (*i).second);
}
}
static void handleName(const TGeoNode* n) {
TGeoVolume* v = n->GetVolume();
TGeoMedium* m = v->GetMedium();
TGeoShape* s = v->GetShape();
string nam;
cout << "Node:'" << n->GetName() << "' Vol:'" << v->GetName() << "' Shape:'" << s->GetName() << "' Medium:'" << m->GetName() << "'" << endl;
}
void Geant4Converter::create(DetElement top) {
G4GeometryInfo geo;
m_dataPtr = &geo;
m_data->clear();
collect(top,geo);
// We do not have to handle defines etc.
// All positions and the like are not really named.
// Hence, start creating the G4 objects for materials, solids and log volumes.
handle(this, geo.materials, &Geant4Converter::handleMaterial);
handle(this, geo.solids, &Geant4Converter::handleSolid);
handle(this, geo.volumes, &Geant4Converter::handleVolume);
// Now place all this stuff appropriately
for(Data::const_reverse_iterator i=m_data->rbegin(); i != m_data->rend(); ++i) {
const Data::mapped_type& v = (*i).second;
for(Data::mapped_type::const_iterator j=v.begin(); j != v.end(); ++j) {
handlePlacement((*j)->GetName(),*j);
handleName(*j);
}
}
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment