diff --git a/DDG4/CMakeLists.txt b/DDG4/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..e5e250834b7ff72971d2be65aeaf9fe5fff663b9 --- /dev/null +++ b/DDG4/CMakeLists.txt @@ -0,0 +1,18 @@ +#---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}) diff --git a/DDG4/include/Geant4Converter.h b/DDG4/include/Geant4Converter.h new file mode 100644 index 0000000000000000000000000000000000000000..4279a51e9e05c1619d64527a608154bebf915a8a --- /dev/null +++ b/DDG4/include/Geant4Converter.h @@ -0,0 +1,125 @@ +// $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 + diff --git a/DDG4/src/Geant4Converter.cpp b/DDG4/src/Geant4Converter.cpp new file mode 100644 index 0000000000000000000000000000000000000000..ada85eb5d99b384c70ca0884d2dabf36b6a1e487 --- /dev/null +++ b/DDG4/src/Geant4Converter.cpp @@ -0,0 +1,313 @@ +// $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); + } + } +}