Geant4Converter.cpp

// $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);
    }
  }
}