TG4MCGeometry.cxx

Go to the documentation of this file.

//------------------------------------------------
// The Geant4 Virtual Monte Carlo package
// Copyright (C) 2007 - 2014 Ivana Hrivnacova
// All rights reserved.
//
// For the licensing terms see geant4_vmc/LICENSE.
// Contact: root-vmc@cern.ch
//-------------------------------------------------
 
 
#include "TG4MCGeometry.h"
#include "TG4G3ControlVector.h"
#include "TG4G3CutVector.h"
#include "TG4G3Units.h"
#include "TG4GeometryServices.h"
#include "TG4Globals.h"
#include "TG4Limits.h"
#include "TG4Medium.h"
 
#ifdef USE_G3TOG4
#include <G3MatTable.hh>
#include <G3toG4.hh>
#endif
 
#include <G4Box.hh>
#include <G4Cons.hh>
#include <G4EllipticalTube.hh>
#include <G4Hype.hh>
#include <G4LogicalBorderSurface.hh>
#include <G4LogicalSkinSurface.hh>
#include <G4LogicalVolume.hh>
#include <G4Material.hh>
#include <G4MaterialPropertiesTable.hh>
#include <G4OpticalSurface.hh>
#include <G4Para.hh>
#include <G4Polycone.hh>
#include <G4Polyhedra.hh>
#include <G4ReflectedSolid.hh>
#include <G4Sphere.hh>
#include <G4SystemOfUnits.hh>
#include <G4Torus.hh>
#include <G4Transform3D.hh>
#include <G4Trap.hh>
#include <G4Trd.hh>
#include <G4Tubs.hh>
#include <G4TwistedTrap.hh>
 
#include <Riostream.h>
#include <TArrayD.h>
#include <TString.h>
 
#ifdef USE_G3TOG4
void G3CLRead(G4String&, char*);
#endif
 
//_____________________________________________________________________________
TG4MCGeometry::TG4MCGeometry()
  : TG4Verbose("g4MCGeometry"), fGeometryServices(0), fMaterialNameVector()
{
 
  fGeometryServices = TG4GeometryServices::Instance();
  if (!fGeometryServices) {
    TG4Globals::Exception("TG4MCGeometry", "TG4MCGeometry",
      "TG4GeometryServices have to be defined first");
  }
}
 
//_____________________________________________________________________________
TG4MCGeometry::~TG4MCGeometry()
{
}
 
//
// public methods
//
 
//_____________________________________________________________________________
void TG4MCGeometry::Material(Int_t& kmat, const char* name, Double_t a,
  Double_t z, Double_t dens, Double_t radl, Double_t absl, Float_t* buf,
  Int_t nwbuf)
{
 
  G4double* bufin = fGeometryServices->CreateG4doubleArray(buf, nwbuf);
  Material(kmat, name, a, z, dens, radl, absl, bufin, nwbuf);
  delete[] bufin;
}
 
#ifdef USE_G3TOG4
//_____________________________________________________________________________
void TG4MCGeometry::Material(Int_t& kmat, const char* name, Double_t a,
  Double_t z, Double_t dens, Double_t radl, Double_t /*absl*/, Double_t* buf,
  Int_t nwbuf)
{
 
  G4String namein = fGeometryServices->CutMaterialName(name);
 
  // create new material only if it does not yet exist
  G4Material* material = fGeometryServices->FindMaterial(a, z, dens);
  if (material) {
    // verbose
    if (VerboseLevel() > 1) {
      G4cout << "!!! Material " << namein << " already exists as "
             << material->GetName() << G4endl;
    }
    G3Mat.put(kmat, material);
  }
  else
    G4gsmate(kmat, namein, a, z, dens, radl, nwbuf, buf);
 
  // save the original material name
  fMaterialNameVector.push_back(namein);
 
  if (nwbuf > 0) {
    TG4Globals::Warning("TG4MCGeometry", "Material",
      "User defined parameters for material " + TString(namein) +
        " are ignored by Geant4.");
  }
}
#else
//_____________________________________________________________________________
void TG4MCGeometry::Material(Int_t&, const char*, Double_t, Double_t, Double_t,
  Double_t, Double_t, Double_t*, Int_t)
{
 
  TG4Globals::Exception("TG4MCGeometry", "Material",
    "This method requires Geant4 installation with G3toG4.");
}
#endif
 
//_____________________________________________________________________________
void TG4MCGeometry::Mixture(Int_t& kmat, const char* name, Float_t* a,
  Float_t* z, Double_t dens, Int_t nlmat, Float_t* wmat)
{
 
  G4double* ain = fGeometryServices->CreateG4doubleArray(a, abs(nlmat));
  G4double* zin = fGeometryServices->CreateG4doubleArray(z, abs(nlmat));
  G4double* wmatin = fGeometryServices->CreateG4doubleArray(wmat, abs(nlmat));
 
  Mixture(kmat, name, ain, zin, dens, nlmat, wmatin);
 
  // !!! in Geant3:
  // After a call with ratios by number (negative number of elements),
  // the ratio array is changed to the ratio by weight, so all successive
  // calls with the same array must specify the number of elements as
  // positive
 
  // wmatin may be modified
  for (G4int i = 0; i < abs(nlmat); i++) wmat[i] = wmatin[i];
 
  delete[] ain;
  delete[] zin;
  delete[] wmatin;
}
 
#ifdef USE_G3TOG4
//_____________________________________________________________________________
void TG4MCGeometry::Mixture(Int_t& kmat, const char* name, Double_t* a,
  Double_t* z, Double_t dens, Int_t nlmat, Double_t* wmat)
{
 
  G4String namein = fGeometryServices->CutMaterialName(name);
 
  // create new material only if it does not yet exist
  G4Material* material =
    fGeometryServices->FindMaterial(a, z, dens, nlmat, wmat);
  if (material) {
    // verbose
    if (VerboseLevel() > 1) {
      G4cout << "!!! Material " << namein << " already exists as "
             << material->GetName() << G4endl;
    }
    G3Mat.put(kmat, material);
  }
  else
    G4gsmixt(kmat, namein, a, z, dens, nlmat, wmat);
 
  // save the original material name
  fMaterialNameVector.push_back(namein);
}
#else
//_____________________________________________________________________________
void TG4MCGeometry::Mixture(
  Int_t&, const char*, Double_t*, Double_t*, Double_t, Int_t, Double_t*)
{
 
  TG4Globals::Exception("TG4MCGeometry", "Mixture",
    "This method requires Geant4 installation with G3toG4.");
}
#endif
 
//_____________________________________________________________________________
void TG4MCGeometry::Medium(Int_t& kmed, const char* name, Int_t nmat,
  Int_t isvol, Int_t ifield, Double_t fieldm, Double_t tmaxfd, Double_t stemax,
  Double_t deemax, Double_t epsil, Double_t stmin, Float_t* ubuf, Int_t nbuf)
{
 
  G4double* bufin = fGeometryServices->CreateG4doubleArray(ubuf, nbuf);
  Medium(kmed, name, nmat, isvol, ifield, fieldm, tmaxfd, stemax, deemax, epsil,
    stmin, bufin, nbuf);
  delete[] bufin;
}
 
#ifdef USE_G3TOG4
//_____________________________________________________________________________
void TG4MCGeometry::Medium(Int_t& kmed, const char* name, Int_t nmat,
  Int_t isvol, Int_t ifield, Double_t fieldm, Double_t tmaxfd, Double_t stemax,
  Double_t deemax, Double_t epsil, Double_t stmin, Double_t* /*ubuf*/,
  Int_t nbuf)
{
 
  G4String namein = fGeometryServices->CutMaterialName(name);
 
  G4gstmed(kmed, name, nmat, isvol, ifield, fieldm, tmaxfd, stemax, deemax,
    epsil, stmin, 0, stemax > 0.);
  // instead of the nbuf argument the bool is passed
  // in order to pass stemax into G4UserLimits
 
  // Add medium to the map (to keep its name which is not stored in G3toG4)
  TG4Medium* medium = fGeometryServices->GetMediumMap()->AddMedium(kmed);
  medium->SetName(name);
  medium->SetIfield(ifield);
 
  if (nbuf > 0) {
    TG4Globals::Warning("TG4MCGeometry", "Medium",
      "User defined parameters for medium " + TString(namein) +
        " are ignored by Geant4.");
  }
}
#else
//_____________________________________________________________________________
void TG4MCGeometry::Medium(Int_t&, const char*, Int_t, Int_t, Int_t, Double_t,
  Double_t, Double_t, Double_t, Double_t, Double_t, Double_t*, Int_t)
{
 
  TG4Globals::Exception("TG4MCGeometry", "Medium",
    "This method requires Geant4 installation with G3toG4.");
}
#endif
 
#ifdef USE_G3TOG4
//_____________________________________________________________________________
void TG4MCGeometry::Matrix(Int_t& krot, Double_t thetaX, Double_t phiX,
  Double_t thetaY, Double_t phiY, Double_t thetaZ, Double_t phiZ)
{
 
  G4gsrotm(krot, thetaX, phiX, thetaY, phiY, thetaZ, phiZ);
}
#else
//_____________________________________________________________________________
void TG4MCGeometry::Matrix(
  Int_t&, Double_t, Double_t, Double_t, Double_t, Double_t, Double_t)
{
 
  TG4Globals::Exception("TG4MCGeometry", "Matrix",
    "This method requires Geant4 installation with G3toG4.");
}
#endif
 
//_____________________________________________________________________________
void TG4MCGeometry::Ggclos()
{
 
#ifdef USE_G3TOG4
  G4ggclos();
#else
  TG4Globals::Exception("TG4MCGeometry", "Ggclos()",
    "This method requires Geant4 installation with G3toG4.");
#endif
}
 
#ifdef USE_G3TOG4
//_____________________________________________________________________________
Int_t TG4MCGeometry::Gsvolu(
  const char* name, const char* shape, Int_t nmed, Double_t* upar, Int_t npar)
{
  // ---
 
  G4gsvolu(fGeometryServices->CutName(name), fGeometryServices->CutName(shape),
    nmed, upar, npar);
  return 0;
}
#else
//_____________________________________________________________________________
Int_t TG4MCGeometry::Gsvolu(const char*, const char*, Int_t, Double_t*, Int_t)
{
 
  TG4Globals::Exception("TG4MCGeometry", "Gsvolu",
    "This method requires Geant4 installation with G3toG4.");
  return 0;
}
#endif
 
//_____________________________________________________________________________
Int_t TG4MCGeometry::Gsvolu(
  const char* name, const char* shape, Int_t nmed, Float_t* upar, Int_t npar)
{
 
  G4double* parin = fGeometryServices->CreateG4doubleArray(upar, npar);
 
  G4int result = Gsvolu(name, shape, nmed, parin, npar);
 
  delete[] parin;
 
  return result;
}
 
#ifdef USE_G3TOG4
//_____________________________________________________________________________
void TG4MCGeometry::Gsdvn(
  const char* name, const char* mother, Int_t ndiv, Int_t iaxis)
{
  //  ---
 
  G4gsdvn(fGeometryServices->CutName(name), fGeometryServices->CutName(mother),
    ndiv, iaxis);
}
#else
//_____________________________________________________________________________
void TG4MCGeometry::Gsdvn(const char*, const char*, Int_t, Int_t)
{
  //  ---
 
  TG4Globals::Exception("TG4MCGeometry", "Gsdvn",
    "This method requires Geant4 installation with G3toG4.");
}
#endif
 
#ifdef USE_G3TOG4
//_____________________________________________________________________________
void TG4MCGeometry::Gsdvn2(const char* name, const char* mother, Int_t ndiv,
  Int_t iaxis, Double_t c0i, Int_t numed)
{
 
  G4gsdvn2(fGeometryServices->CutName(name), fGeometryServices->CutName(mother),
    ndiv, iaxis, c0i, numed);
}
#else
//_____________________________________________________________________________
void TG4MCGeometry::Gsdvn2(
  const char*, const char*, Int_t, Int_t, Double_t, Int_t)
{
 
  TG4Globals::Exception("TG4MCGeometry", "Gsdvn2",
    "This method requires Geant4 installation with G3toG4.");
}
#endif
 
#ifdef USE_G3TOG4
//_____________________________________________________________________________
void TG4MCGeometry::Gsdvt(const char* name, const char* mother, Double_t step,
  Int_t iaxis, Int_t numed, Int_t ndvmx)
{
 
  G4gsdvt(fGeometryServices->CutName(name), fGeometryServices->CutName(mother),
    step, iaxis, numed, ndvmx);
}
#else
//_____________________________________________________________________________
void TG4MCGeometry::Gsdvt(
  const char*, const char*, Double_t, Int_t, Int_t, Int_t)
{
 
  TG4Globals::Exception("TG4MCGeometry", "Gsdvt",
    "This method requires Geant4 installation with G3toG4.");
}
#endif
 
#ifdef USE_G3TOG4
//_____________________________________________________________________________
void TG4MCGeometry::Gsdvt2(const char* name, const char* mother, Double_t step,
  Int_t iaxis, Double_t c0, Int_t numed, Int_t ndvmx)
{
 
  G4gsdvt2(fGeometryServices->CutName(name), fGeometryServices->CutName(mother),
    step, iaxis, c0, numed, ndvmx);
}
#else
//_____________________________________________________________________________
void TG4MCGeometry::Gsdvt2(
  const char*, const char*, Double_t, Int_t, Double_t, Int_t, Int_t)
{
 
  TG4Globals::Exception("TG4MCGeometry", "Gsdvt",
    "This method requires Geant4 installation with G3toG4.");
}
#endif
 
//_____________________________________________________________________________
void TG4MCGeometry::Gsord(const char* /*name*/, Int_t /*iax*/)
{
 
  TG4Globals::Warning("TG4MCGeometry", "Gsord", "Dummy method.");
}
 
#ifdef USE_G3TOG4
//_____________________________________________________________________________
void TG4MCGeometry::Gspos(const char* vname, Int_t num, const char* vmoth,
  Double_t x, Double_t y, Double_t z, Int_t irot, const char* vonly)
{
  // ---
 
  G4gspos(fGeometryServices->CutName(vname), ++num,
    fGeometryServices->CutName(vmoth), x, y, z, irot, vonly);
}
#else
//_____________________________________________________________________________
void TG4MCGeometry::Gspos(const char*, Int_t, const char*, Double_t, Double_t,
  Double_t, Int_t, const char*)
{
  // ---
 
  TG4Globals::Exception("TG4MCGeometry", "Gspos",
    "This method requires Geant4 installation with G3toG4.");
}
#endif
 
#ifdef USE_G3TOG4
//_____________________________________________________________________________
void TG4MCGeometry::Gsposp(const char* name, Int_t nr, const char* mother,
  Double_t x, Double_t y, Double_t z, Int_t irot, const char* konly,
  Double_t* upar, Int_t np)
{
 
  G4gsposp(fGeometryServices->CutName(name), ++nr,
    fGeometryServices->CutName(mother), x, y, z, irot, konly, upar, np);
}
#else
//_____________________________________________________________________________
void TG4MCGeometry::Gsposp(const char*, Int_t, const char*, Double_t, Double_t,
  Double_t, Int_t, const char*, Double_t*, Int_t)
{
 
  TG4Globals::Exception("TG4MCGeometry", "Gsposp",
    "This method requires Geant4 installation with G3toG4.");
}
#endif
 
//_____________________________________________________________________________
void TG4MCGeometry::Gsposp(const char* name, Int_t nr, const char* mother,
  Double_t x, Double_t y, Double_t z, Int_t irot, const char* konly,
  Float_t* upar, Int_t np)
{
 
  G4double* parin = fGeometryServices->CreateG4doubleArray(upar, np);
  Gsposp(name, nr, mother, x, y, z, irot, konly, parin, np);
  delete[] parin;
}
 
#ifdef USE_G3TOG4
//_____________________________________________________________________________
void TG4MCGeometry::Gsbool(const char* onlyVolName, const char* manyVolName)
{
  // ---
 
  G4gsbool(onlyVolName, manyVolName);
}
#else
//_____________________________________________________________________________
void TG4MCGeometry::Gsbool(const char*, const char*)
{
  // ---
 
  TG4Globals::Exception("TG4MCGeometry", "Gsbool",
    "This method requires Geant4 installation with G3toG4.");
}
#endif
 
//_____________________________________________________________________________
Bool_t TG4MCGeometry::GetTransformation(
  const TString& volumePath, TGeoHMatrix& matrix)
{
 
  G4String path = volumePath.Data();
 
  G4String volName;
  G4int copyNo;
  path = fGeometryServices->CutVolumePath(path, volName, copyNo);
 
  // Get the first volume
  G4VPhysicalVolume* pvTop = 0;
  if (fGeometryServices->GetWorld()->GetName() == volName) {
    // Check world volume first
    pvTop = fGeometryServices->GetWorld();
  }
  else
    pvTop = fGeometryServices->FindPhysicalVolume(volName, copyNo, true);
 
  if (!pvTop) {
    TG4Globals::Warning("TG4MCGeometry",
      "GetTransformation:", "Top volume in " + volumePath + " does not exist.");
    return false;
  }
 
  // Calculate transformation
  //
  G4VPhysicalVolume* pvMother = pvTop;
  G4Transform3D transform;
 
  while (path.length() > 0) {
    // Extract next volume name & copyNo
    path = fGeometryServices->CutVolumePath(path, volName, copyNo);
 
    // Find daughter
    G4VPhysicalVolume* pvDaughter = fGeometryServices->FindDaughter(
      volName, copyNo, pvMother->GetLogicalVolume(), true);
 
    if (!pvDaughter) {
      TG4Globals::Warning("TG4MCGeometry", "GetTransformation",
        "Daughter volume " + TString(volName) + " in " + volumePath +
          " does not exist.");
      return false;
    }
 
    transform = transform * G4Transform3D(*pvDaughter->GetObjectRotation(),
                              pvDaughter->GetObjectTranslation());
    pvMother = pvDaughter;
  }
 
  // put transform in TGeoHMatrix here
  fGeometryServices->Convert(transform, matrix);
  return true;
}
 
//_____________________________________________________________________________
Bool_t TG4MCGeometry::GetShape(
  const TString& volumePath, TString& shapeType, TArrayD& par)
{
 
  // Get volume & copyNo
  G4String path = volumePath.Data();
  G4int last1 = path.rfind('/');
  G4int last2 = path.rfind('_');
  G4String volName = path.substr(last1 + 1, last2 - last1 - 1);
  G4String copyNoStr = path.substr(last2 + 1, path.length() - last2);
  std::istringstream in(copyNoStr);
  G4int copyNo;
  in >> copyNo;
 
  // Find physical volume
  G4VPhysicalVolume* pv =
    fGeometryServices->FindPhysicalVolume(volName, copyNo, true);
 
  if (!pv) {
    TG4Globals::Warning("TG4MCGeometry", "GetShape",
      "Physical volume " + volumePath + " does not exist.");
    return false;
  }
 
  // Get solid
  G4VSolid* solid = pv->GetLogicalVolume()->GetSolid();
 
  // Get constituent solid if reflected solid
  Bool_t isReflected = false;
  G4ReflectedSolid* reflSolid = dynamic_cast<G4ReflectedSolid*>(solid);
  if (reflSolid) {
    solid = reflSolid->GetConstituentMovedSolid();
    isReflected = true;
  }
 
  Int_t npar;
 
  if (solid->GetEntityType() == "G4Box") {
    shapeType = "BOX";
    npar = 3;
    par.Set(npar);
    G4Box* box = (G4Box*)solid;
    par.AddAt(box->GetXHalfLength() / cm, 0);
    par.AddAt(box->GetYHalfLength() / cm, 1);
    par.AddAt(box->GetZHalfLength() / cm, 2);
    return kTRUE;
  }
 
  if (solid->GetEntityType() == "G4Cons") {
    shapeType = "CONS";
    npar = 7;
    par.Set(npar);
    G4Cons* cons = (G4Cons*)solid;
    par.AddAt(cons->GetZHalfLength() / cm, 0);
    if (!isReflected) {
      par.AddAt(cons->GetInnerRadiusMinusZ() / cm, 1);
      par.AddAt(cons->GetOuterRadiusMinusZ() / cm, 2);
      par.AddAt(cons->GetInnerRadiusPlusZ() / cm, 3);
      par.AddAt(cons->GetOuterRadiusPlusZ() / cm, 4);
    }
    else {
      par.AddAt(cons->GetInnerRadiusPlusZ() / cm, 1);
      par.AddAt(cons->GetOuterRadiusPlusZ() / cm, 2);
      par.AddAt(cons->GetInnerRadiusMinusZ() / cm, 3);
      par.AddAt(cons->GetOuterRadiusMinusZ() / cm, 4);
    }
    par.AddAt(cons->GetStartPhiAngle() / deg, 5);
    par.AddAt(cons->GetDeltaPhiAngle() / deg, 6);
    return kTRUE;
  }
 
  if (solid->GetEntityType() == "G4EllipticalTube") {
    shapeType = "ELTU";
    npar = 3;
    par.Set(npar);
    G4EllipticalTube* eltu = (G4EllipticalTube*)solid;
    par.AddAt(eltu->GetDx() / cm, 0);
    par.AddAt(eltu->GetDy() / cm, 1);
    par.AddAt(eltu->GetDz() / cm, 2);
    return kTRUE;
  }
 
  // Add to G3toG4, VGM
  if (solid->GetEntityType() == "G4Hype") {
    shapeType = "HYPE";
    npar = 5;
    par.Set(npar);
    G4Hype* hype = (G4Hype*)solid;
    par.AddAt(hype->GetInnerRadius() / cm, 0);
    par.AddAt(hype->GetOuterRadius() / cm, 1);
    par.AddAt(hype->GetZHalfLength() / cm, 2);
    par.AddAt(hype->GetInnerStereo(), 3); // ?? unit
    par.AddAt(hype->GetOuterStereo(), 4); // ?? unit
    return kTRUE;
  }
 
  // Add G4Orb
 
  if (solid->GetEntityType() == "G4Para") {
    shapeType = "PARA";
    npar = 6;
    par.Set(npar);
    G4Para* para = (G4Para*)solid;
    par.AddAt(para->GetXHalfLength() / cm, 0);
    par.AddAt(para->GetYHalfLength() / cm, 1);
    par.AddAt(para->GetZHalfLength() / cm, 2);
    par.AddAt(atan(para->GetTanAlpha()) / deg, 3);
    if (!isReflected)
      par.AddAt(para->GetSymAxis().theta() / deg, 4);
    else
      par.AddAt((M_PI - para->GetSymAxis().theta()) / deg, 4);
    par.AddAt(para->GetSymAxis().phi() / deg, 5);
    return kTRUE;
  }
 
  if (solid->GetEntityType() == "G4Polycone") {
    shapeType = "PCON";
    G4Polycone* pcon = (G4Polycone*)solid;
    Int_t nz = pcon->GetOriginalParameters()->Num_z_planes;
    const Double_t* rmin = pcon->GetOriginalParameters()->Rmin;
    const Double_t* rmax = pcon->GetOriginalParameters()->Rmax;
    const Double_t* z = pcon->GetOriginalParameters()->Z_values;
    npar = 3 + 3 * nz;
    par.Set(npar);
    par.AddAt(pcon->GetStartPhi() / deg, 0);
    par.AddAt((pcon->GetEndPhi() - pcon->GetStartPhi()) / deg, 1);
    par.AddAt(pcon->GetOriginalParameters()->Num_z_planes, 2);
    for (Int_t i = 0; i < nz; i++) {
      if (!isReflected)
        par.AddAt(z[i] / cm, 3 + 3 * i);
      else
        par.AddAt(-z[i] / cm, 3 + 3 * i);
      par.AddAt(rmin[i] / cm, 3 + 3 * i + 1);
      par.AddAt(rmax[i] / cm, 3 + 3 * i + 2);
    }
    return kTRUE;
  }
 
  if (solid->GetEntityType() == "G4Polyhedra") {
    shapeType = "PGON";
    G4Polyhedra* pgon = (G4Polyhedra*)solid;
    Int_t nz = pgon->GetOriginalParameters()->Num_z_planes;
    const Double_t* rmin = pgon->GetOriginalParameters()->Rmin;
    const Double_t* rmax = pgon->GetOriginalParameters()->Rmax;
    const Double_t* z = pgon->GetOriginalParameters()->Z_values;
    npar = 4 + 3 * nz;
    par.Set(npar);
    par.AddAt(pgon->GetStartPhi() / deg, 0);
    par.AddAt((pgon->GetEndPhi() - pgon->GetStartPhi()) / deg, 1);
    par.AddAt(pgon->GetOriginalParameters()->numSide, 2);
    par.AddAt(pgon->GetOriginalParameters()->Num_z_planes, 3);
    for (Int_t i = 0; i < nz; i++) {
      if (!isReflected)
        par.AddAt(z[i] / cm, 4 + 3 * i);
      else
        par.AddAt(-z[i] / cm, 4 + 3 * i);
      par.AddAt(rmin[i] / cm, 4 + 3 * i + 1);
      par.AddAt(rmax[i] / cm, 4 + 3 * i + 2);
    }
    return kTRUE;
  }
 
  if (solid->GetEntityType() == "G4Sphere") {
    shapeType = "SPHE";
    npar = 6;
    par.Set(npar);
    G4Sphere* sphe = (G4Sphere*)solid;
    par.AddAt(sphe->GetInnerRadius() / cm, 0);
    par.AddAt(sphe->GetOuterRadius() / cm, 1);
    if (!isReflected)
      par.AddAt(sphe->GetStartThetaAngle() / deg, 2);
    else
      par.AddAt((M_PI - sphe->GetStartThetaAngle()) / deg, 2);
    par.AddAt(sphe->GetStartPhiAngle() / deg, 3);
    par.AddAt(sphe->GetStartPhiAngle() / deg, 4);
    par.AddAt(sphe->GetDeltaPhiAngle() / deg, 5);
    return kTRUE;
  }
 
  // Add G4Torus
 
  if (solid->GetEntityType() == "G4Trap") {
    shapeType = "TRAP";
    npar = 11;
    par.Set(npar);
    G4Trap* trap = (G4Trap*)solid;
    par.AddAt(trap->GetZHalfLength() / cm, 0);
    if (!isReflected) {
      par.AddAt(trap->GetSymAxis().theta() / deg, 1);
      par.AddAt(trap->GetSymAxis().phi() / deg, 2);
      par.AddAt(trap->GetYHalfLength1() / cm, 3);
      par.AddAt(trap->GetXHalfLength1() / cm, 4);
      par.AddAt(trap->GetXHalfLength2() / cm, 5);
      par.AddAt(atan(trap->GetTanAlpha1()) / deg, 6);
      par.AddAt(trap->GetYHalfLength2() / cm, 7);
      par.AddAt(trap->GetXHalfLength3() / cm, 8);
      par.AddAt(trap->GetXHalfLength4() / cm, 9);
    }
    else {
      par.AddAt(-(trap->GetSymAxis().theta()) / deg, 1);
      par.AddAt(-(trap->GetSymAxis().phi()) / deg, 2);
      par.AddAt(trap->GetYHalfLength2() / cm, 3);
      par.AddAt(trap->GetXHalfLength3() / cm, 4);
      par.AddAt(trap->GetXHalfLength4() / cm, 5);
      par.AddAt(atan(trap->GetTanAlpha2()) / deg, 6);
      par.AddAt(trap->GetYHalfLength1() / cm, 7);
      par.AddAt(trap->GetXHalfLength1() / cm, 8);
      par.AddAt(trap->GetXHalfLength2() / cm, 9);
      par.AddAt(atan(trap->GetTanAlpha1()) / deg, 10);
    }
    return kTRUE;
  }
 
  if (solid->GetEntityType() == "G4Trd") {
    shapeType = "TRD2";
    npar = 5;
    par.Set(npar);
    G4Trd* trd2 = (G4Trd*)solid;
    if (!isReflected) {
      par.AddAt(trd2->GetXHalfLength1() / cm, 0);
      par.AddAt(trd2->GetXHalfLength2() / cm, 1);
      par.AddAt(trd2->GetYHalfLength1() / cm, 2);
      par.AddAt(trd2->GetYHalfLength2() / cm, 3);
    }
    else {
      par.AddAt(trd2->GetXHalfLength2() / cm, 0);
      par.AddAt(trd2->GetXHalfLength1() / cm, 1);
      par.AddAt(trd2->GetYHalfLength2() / cm, 2);
      par.AddAt(trd2->GetYHalfLength1() / cm, 3);
    }
    par.AddAt(trd2->GetZHalfLength() / cm, 4);
    return kTRUE;
  }
 
  if (solid->GetEntityType() == "G4Tubs") {
    shapeType = "TUBS";
    npar = 5;
    par.Set(npar);
    G4Tubs* tubs = (G4Tubs*)solid;
    par.AddAt(tubs->GetInnerRadius() / cm, 0);
    par.AddAt(tubs->GetOuterRadius() / cm, 1);
    par.AddAt(tubs->GetZHalfLength() / cm, 2);
    par.AddAt(tubs->GetStartPhiAngle() / deg, 3);
    par.AddAt(tubs->GetDeltaPhiAngle() / deg, 4);
    return kTRUE;
  }
 
  // To do: add to VGM
  if (solid->GetEntityType() == "G4TwistedTrap") {
    shapeType = "GTRA";
    npar = 12;
    par.Set(npar);
    G4TwistedTrap* trap = (G4TwistedTrap*)solid;
    par.AddAt(trap->GetZHalfLength() / cm, 0);
    if (!isReflected) {
      par.AddAt(trap->GetPolarAngleTheta() / deg, 1);
      par.AddAt(trap->GetAzimuthalAnglePhi() / deg, 2);
      par.AddAt(trap->GetY1HalfLength() / cm, 3);
      par.AddAt(trap->GetX1HalfLength() / cm, 4);
      par.AddAt(trap->GetX2HalfLength() / cm, 5);
      par.AddAt(trap->GetTiltAngleAlpha() / deg, 6);
      par.AddAt(trap->GetY2HalfLength() / cm, 7);
      par.AddAt(trap->GetX3HalfLength() / cm, 8);
      par.AddAt(trap->GetX4HalfLength() / cm, 9);
      par.AddAt(trap->GetTiltAngleAlpha() / deg, 10);
    }
    else {
      par.AddAt(-trap->GetPolarAngleTheta() / deg, 1);
      par.AddAt(-trap->GetAzimuthalAnglePhi() / deg, 2);
      par.AddAt(trap->GetY2HalfLength() / cm, 3);
      par.AddAt(trap->GetX3HalfLength() / cm, 4);
      par.AddAt(trap->GetX4HalfLength() / cm, 5);
      par.AddAt(trap->GetTiltAngleAlpha() / deg, 6);
      par.AddAt(trap->GetY1HalfLength() / cm, 7);
      par.AddAt(trap->GetX1HalfLength() / cm, 8);
      par.AddAt(trap->GetX2HalfLength() / cm, 9);
      par.AddAt(trap->GetTiltAngleAlpha() / deg, 10);
    }
    par.AddAt(trap->GetPhiTwist() / deg, 11);
    return kTRUE;
  }
 
  TG4Globals::Warning("TG4MCGeometry", "GetShape",
    "Shape " + TString(solid->GetEntityType()) + " not implemented.");
  return false;
}
 
//_____________________________________________________________________________
Bool_t TG4MCGeometry::GetMaterial(const TString& volumeName, TString& name,
  Int_t& imat, Double_t& a, Double_t& z, Double_t& density, Double_t& radl,
  Double_t& inter, TArrayD& par)
{
 
  // Get volume
  G4LogicalVolume* lv =
    fGeometryServices->FindLogicalVolume(volumeName.Data(), true);
  if (!lv) {
    TG4Globals::Warning("TG4MCGeometry", "GetMaterial",
      "Logical volume " + volumeName + " does not exist.");
    return false;
  }
 
  G4Material* material = lv->GetMaterial();
  imat = material->GetIndex();
  name = material->GetName();
  a = fGeometryServices->GetEffA(material);
  z = fGeometryServices->GetEffZ(material);
 
  density = material->GetDensity();
  density /= TG4G3Units::MassDensity();
 
  radl = material->GetRadlen();
  radl /= TG4G3Units::Length();
 
  // the following parameters are not defined in Geant4
  inter = 0.;
  par.Set(0);
  return true;
}
 
//_____________________________________________________________________________
Bool_t TG4MCGeometry::GetMedium(const TString& volumeName, TString& name,
  Int_t& imed, Int_t& nmat, Int_t& isvol, Int_t& ifield, Double_t& fieldm,
  Double_t& tmaxfd, Double_t& stemax, Double_t& deemax, Double_t& epsil,
  Double_t& stmin, TArrayD& par)
{
 
  // Get volume
  G4LogicalVolume* lv =
    fGeometryServices->FindLogicalVolume(volumeName.Data(), true);
  if (!lv) {
    TG4Globals::Warning("TG4MCGeometry", "GetMaterial",
      "Logical volume " + volumeName + " does not exist.");
    return false;
  }
 
  imed = fGeometryServices->GetMediumId(lv);
  nmat = lv->GetMaterial()->GetIndex();
 
  TG4Limits* limits = fGeometryServices->GetLimits(lv->GetUserLimits());
  if (limits) {
    name = limits->GetName();
    stemax = limits->GetMaxUserStep();
  }
  else
    stemax = 0.;
 
  // the following parameters are not defined in Geant4
  isvol = 0;
  ifield = 0;
  fieldm = 0;
  tmaxfd = 0.;
  deemax = 0.;
  epsil = 0.;
  stmin = 0.;
  par.Set(0);
 
  return true;
}
 
//_____________________________________________________________________________
Int_t TG4MCGeometry::MediumId(const Text_t* mediumName) const
{
 
  TG4MediumMap* mediumMap = TG4GeometryServices::Instance()->GetMediumMap();
  TG4Medium* medium = mediumMap->GetMedium(G4String(mediumName), false);
 
  if (!medium) {
    TG4Globals::Warning("TG4MCGeometry", "MediumId",
      "Medium " + TString(mediumName) + " not found.");
    return 0;
  }
 
  return medium->GetID();
}
 
//
// Not implemented functions
//
 
//_____________________________________________________________________________
Int_t TG4MCGeometry::VolId(const Text_t* /*volName*/) const
{
 
  TG4Globals::Exception("TG4MCGeometry", "VolId", "Not implemented.");
  return 0;
}
 
//_____________________________________________________________________________
const char* TG4MCGeometry::VolName(Int_t /*id*/) const
{
 
  TG4Globals::Exception("TG4MCGeometry", "VolName", "Not implemented.");
  return "";
}
 
//_____________________________________________________________________________
Int_t TG4MCGeometry::NofVolumes() const
{
 
  TG4Globals::Exception("TG4MCGeometry", "NofVolumes", "Not implemented.");
  return 0;
}
 
//_____________________________________________________________________________
Int_t TG4MCGeometry::NofVolDaughters(const char* /*volName*/) const
{
 
  TG4Globals::Exception("TG4MCGeometry", "NofVolDaughters", "Not implemented.");
  return 0;
}
 
//_____________________________________________________________________________
const char* TG4MCGeometry::VolDaughterName(
  const char* /*volName*/, Int_t /*i*/) const
{
 
  TG4Globals::Exception("TG4MCGeometry", "VolDaughterName", "Not implemented.");
  return "";
}
 
//_____________________________________________________________________________
Int_t TG4MCGeometry::VolDaughterCopyNo(
  const char* /*volName*/, Int_t /*i*/) const
{
 
  TG4Globals::Exception(
    "TG4MCGeometry", "VolDaughterCopyNo", "Not implemented.");
  return 0;
}
 
//_____________________________________________________________________________
Int_t TG4MCGeometry::VolId2Mate(Int_t /*id*/) const
{
 
  TG4Globals::Exception("TG4MCGeometry", "VolId2Mate", "Not implemented.");
  return 0;
}