Newer
Older
Markus Frank
committed
//==========================================================================
Markus Frank
committed
//--------------------------------------------------------------------------
// Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
Markus Frank
committed
// All rights reserved.
Markus Frank
committed
// For the licensing terms see $DD4hepINSTALL/LICENSE.
// For the list of contributors see $DD4hepINSTALL/doc/CREDITS.
Markus Frank
committed
// Author : M.Frank
//
//==========================================================================
Markus Frank
committed
// Framework include files
#include "DD4hep/Plugins.h"
Markus Frank
committed
#include "DD4hep/Printout.h"
#include "DD4hep/DD4hepUnits.h"
#include "DD4hep/PropertyTable.h"
#include "DD4hep/detail/ObjectsInterna.h"
#include "DD4hep/detail/DetectorInterna.h"
#include "DDG4/Geant4Field.h"
#include "DDG4/Geant4Converter.h"
// ROOT includes
#include "TROOT.h"
#include "TColor.h"
Markus Frank
committed
#include "TGeoNode.h"
#include "TGeoShape.h"
#include "TGeoCone.h"
#include "TGeoEltu.h"
#include "TGeoPgon.h"
#include "TGeoSphere.h"
#include "TGeoTorus.h"
#include "TGeoTrd1.h"
#include "TGeoTrd2.h"
#include "TGeoArb8.h"
#include "TGeoMatrix.h"
#include "TGeoBoolNode.h"
Markus Frank
committed
#include "TGeoParaboloid.h"
Markus Frank
committed
#include "TGeoShapeAssembly.h"
#include "TGeoManager.h"
#include "TClass.h"
#include "TMath.h"
#include "G4VisAttributes.hh"
#include "G4ProductionCuts.hh"
#include "G4VUserRegionInformation.hh"
#include "G4Element.hh"
#include "G4Box.hh"
#include "G4Trd.hh"
#include "G4Polycone.hh"
#include "G4Polyhedra.hh"
#include "G4UnionSolid.hh"
#include "G4Ellipsoid.hh"
#include "G4GenericTrap.hh"
#include "G4ExtrudedSolid.hh"
#include "G4ReflectedSolid.hh"
#include "G4EllipticalTube.hh"
#include "G4SubtractionSolid.hh"
#include "G4IntersectionSolid.hh"
#include "G4Region.hh"
#include "G4UserLimits.hh"
#include "G4LogicalVolume.hh"
#include "G4Material.hh"
#include "G4Element.hh"
#include "G4Isotope.hh"
#include "G4Transform3D.hh"
#include "G4ThreeVector.hh"
#include "G4PVPlacement.hh"
#include "G4ElectroMagneticField.hh"
#include "G4FieldManager.hh"
#include "G4ReflectionFactory.hh"
#include "G4OpticalSurface.hh"
#include "G4LogicalSkinSurface.hh"
#include "G4LogicalBorderSurface.hh"
#include "G4MaterialPropertiesTable.hh"
#include "G4MaterialPropertiesIndex.hh"
#include "CLHEP/Units/SystemOfUnits.h"
Markus Frank
committed
// C/C++ include files
Markus Frank
committed
#include <iostream>
#include <iomanip>
#include <sstream>
namespace units = dd4hep;
using namespace dd4hep::detail;
using namespace dd4hep::sim;
using namespace dd4hep;
#include "DDG4/Geant4AssemblyVolume.h"
#include "DD4hep/DetectorTools.h"
#include "G4RotationMatrix.hh"
#include "G4AffineTransform.hh"
#include "G4LogicalVolume.hh"
#include "G4VPhysicalVolume.hh"
#include "G4ReflectionFactory.hh"
static const double CM_2_MM = (CLHEP::centimeter/dd4hep::centimeter);
void Geant4AssemblyVolume::imprint(Geant4GeometryInfo& info,
Markus Frank
committed
const TGeoNode* parent,
Chain chain,
Geant4AssemblyVolume* pAssembly,
G4LogicalVolume* pMotherLV,
G4Transform3D& transformation,
G4int copyNumBase,
G4bool surfCheck )
TGeoVolume* vol = parent->GetVolume();
static int level=0;
++level;
Markus Frank
committed
unsigned int numberOfDaughters;
if( copyNumBase == 0 )
Markus Frank
committed
{
numberOfDaughters = pMotherLV->GetNoDaughters();
}
Markus Frank
committed
{
numberOfDaughters = copyNumBase;
}
// We start from the first available index
numberOfDaughters++;
ImprintsCountPlus();
Markus Frank
committed
vector<G4AssemblyTriplet> triplets = pAssembly->fTriplets;
//cout << " Assembly:" << detail::tools::placementPath(chain) << endl;
for( unsigned int i = 0; i < triplets.size(); i++ ) {
const TGeoNode* node = pAssembly->m_entries[i];
Chain new_chain = chain;
new_chain.emplace_back(node);
//cout << " Assembly: Entry: " << detail::tools::placementPath(new_chain) << endl;
G4Transform3D Ta( *(triplets[i].GetRotation()),
triplets[i].GetTranslation() );
if ( triplets[i].IsReflection() ) { Ta = Ta * G4ReflectZ3D(); }
G4Transform3D Tfinal = transformation * Ta;
Markus Frank
committed
Markus Frank
committed
{
// Generate the unique name for the next PV instance
// The name has format:
//
// av_WWW_impr_XXX_YYY_ZZZ
// where the fields mean:
// WWW - assembly volume instance number
// XXX - assembly volume imprint number
// YYY - the name of a log. volume we want to make a placement of
// ZZZ - the log. volume index inside the assembly volume
//
stringstream pvName;
Markus Frank
committed
pvName << "av_"
<< GetAssemblyID()
<< "_impr_"
<< GetImprintsCount()
<< "_"
<< triplets[i].GetVolume()->GetName().c_str()
<< "_pv_"
<< i
<< ends;
Markus Frank
committed
// Generate a new physical volume instance inside a mother
// (as we allow 3D transformation use G4ReflectionFactory to
// take into account eventual reflection)
//
G4PhysicalVolumesPair pvPlaced
= G4ReflectionFactory::Instance()->Place( Tfinal,
pvName.str().c_str(),
triplets[i].GetVolume(),
pMotherLV,
false,
numberOfDaughters + i,
surfCheck );
// Register the physical volume created by us so we can delete it later
//
//fPVStore.emplace_back( pvPlaced.first );
info.g4VolumeImprints[vol].emplace_back(new_chain,pvPlaced.first);
#if 0
Markus Frank
committed
cout << " Assembly:Parent:" << parent->GetName() << " " << node->GetName()
<< " " << (void*)node << " G4:" << pvName.str() << " Daughter:"
<< detail::tools::placementPath(new_chain) << endl;
Markus Frank
committed
cout << endl;
#endif
Markus Frank
committed
if ( pvPlaced.second ) {
G4Exception("G4AssemblyVolume::MakeImprint(..)", "GeomVol0003", FatalException,
"Fancy construct popping new mother from the stack!");
//fPVStore.emplace_back( pvPlaced.second );
Markus Frank
committed
}
else if ( triplets[i].GetAssembly() ) {
// Place volumes in this assembly with composed transformation
Markus Frank
committed
imprint(info, parent, new_chain, (Geant4AssemblyVolume*)triplets[i].GetAssembly(), pMotherLV, Tfinal, i*100+copyNumBase, surfCheck );
else {
--level;
G4Exception("G4AssemblyVolume::MakeImprint(..)",
"GeomVol0003", FatalException,
"Triplet has no volume and no assembly");
Markus Frank
committed
}
//cout << "Imprinted assembly level:" << level << " in mother:" << pMotherLV->GetName() << endl;
Markus Frank
committed
}
namespace {
static string indent = "";
static Double_t s_identity_rot[] = { 1., 0., 0., 0., 1., 0., 0., 0., 1. };
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,
Markus Frank
committed
double ZZ, double DZ)
: G4Transform3D(XX, XY, XZ, DX, YX, YY, YZ, DY, ZX, ZY, ZZ, DZ) {
Markus Frank
committed
MyTransform3D(const double* t, const double* r = s_identity_rot)
: G4Transform3D(r[0],r[1],r[2],t[0]*CM_2_MM,r[3],r[4],r[5],t[1]*CM_2_MM,r[6],r[7],r[8],t[2]*CM_2_MM) {
}
Markus Frank
committed
#if 0 // warning: unused function 'handleName' [-Wunused-function]
void handleName(const TGeoNode* n) {
TGeoVolume* v = n->GetVolume();
TGeoMedium* m = v->GetMedium();
printout(DEBUG, "G4", "TGeoNode:'%s' Vol:'%s' Shape:'%s' Medium:'%s'", n->GetName(), v->GetName(), s->GetName(),
Markus Frank
committed
m->GetName());
Markus Frank
committed
#endif
class G4UserRegionInformation : public G4VUserRegionInformation {
public:
Region region;
double threshold;
bool storeSecondaries;
G4UserRegionInformation()
Markus Frank
committed
: threshold(0.0), storeSecondaries(false) {
}
virtual ~G4UserRegionInformation() {
}
if (region.isValid())
printout(DEBUG, "Region", "Name:%s", region.name());
pair<double,double> g4PropertyConversion(int index) {
switch(index) {
case kRINDEX: return make_pair(CLHEP::keV/units::keV, 1.0);
case kREFLECTIVITY: return make_pair(CLHEP::keV/units::keV, 1.0);
case kREALRINDEX: return make_pair(CLHEP::keV/units::keV, 1.0);
case kIMAGINARYRINDEX: return make_pair(CLHEP::keV/units::keV, 1.0);
case kEFFICIENCY: return make_pair(CLHEP::keV/units::keV, 1.0);
case kTRANSMITTANCE: return make_pair(CLHEP::keV/units::keV, 1.0);
case kSPECULARLOBECONSTANT: return make_pair(CLHEP::keV/units::keV, 1.0);
case kSPECULARSPIKECONSTANT: return make_pair(CLHEP::keV/units::keV, 1.0);
case kBACKSCATTERCONSTANT: return make_pair(CLHEP::keV/units::keV, 1.0);
case kGROUPVEL: return make_pair(CLHEP::keV/units::keV, (CLHEP::m/CLHEP::s)/(units::m/units::s)); // meter/second
case kMIEHG: return make_pair(CLHEP::keV/units::keV, CLHEP::m/units::m);
case kRAYLEIGH: return make_pair(CLHEP::keV/units::keV, CLHEP::m/units::m); // ??? says its a length
case kWLSCOMPONENT: return make_pair(CLHEP::keV/units::keV, 1.0);
case kWLSABSLENGTH: return make_pair(CLHEP::keV/units::keV, CLHEP::m/units::m);
case kABSLENGTH: return make_pair(CLHEP::keV/units::keV, CLHEP::m/units::m);
case kFASTCOMPONENT: return make_pair(CLHEP::keV/units::keV, 1.0);
case kSLOWCOMPONENT: return make_pair(CLHEP::keV/units::keV, 1.0);
case kPROTONSCINTILLATIONYIELD: return make_pair(CLHEP::keV/units::keV, units::keV/CLHEP::keV); // Yields: 1/energy
case kDEUTERONSCINTILLATIONYIELD: return make_pair(CLHEP::keV/units::keV, units::keV/CLHEP::keV);
case kTRITONSCINTILLATIONYIELD: return make_pair(CLHEP::keV/units::keV, units::keV/CLHEP::keV);
case kALPHASCINTILLATIONYIELD: return make_pair(CLHEP::keV/units::keV, units::keV/CLHEP::keV);
case kIONSCINTILLATIONYIELD: return make_pair(CLHEP::keV/units::keV, units::keV/CLHEP::keV);
case kELECTRONSCINTILLATIONYIELD: return make_pair(CLHEP::keV/units::keV, units::keV/CLHEP::keV);
default:
break;
}
printout(FATAL,"Geant4Converter", "+++ Cannot convert material property with index: %d", index);
#else
printout(FATAL,"Geant4Converter", "+++ Cannot convert material property with index: %d [Need Geant4 > 10.03]", index);
#endif
return make_pair(0e0,0e0);
}
double g4ConstPropertyConversion(int index) {
case kSURFACEROUGHNESS: return CLHEP::m/units::m; // Length
case kISOTHERMAL_COMPRESSIBILITY: return (CLHEP::m3/CLHEP::keV)/(units::m3/CLHEP::keV); // Volume/Energy
case kRS_SCALE_FACTOR: return 1.0; // ??
case kWLSMEANNUMBERPHOTONS: return 1.0; // ??
case kWLSTIMECONSTANT: return CLHEP::second/units::second; // Time
case kMIEHG_FORWARD: return 1.0;
case kMIEHG_BACKWARD: return 1.0;
case kMIEHG_FORWARD_RATIO: return 1.0;
case kSCINTILLATIONYIELD: return units::keV/CLHEP::keV; // Energy
case kRESOLUTIONSCALE: return 1.0;
case kFASTTIMECONSTANT: return CLHEP::second/units::second; // Time
case kFASTSCINTILLATIONRISETIME: return CLHEP::second/units::second; // Time
case kSLOWTIMECONSTANT: return CLHEP::second/units::second; // Time
case kSLOWSCINTILLATIONRISETIME: return CLHEP::second/units::second; // Time
case kFERMIPOT: return CLHEP::keV/units::keV; // Energy
case kDIFFUSION: return 1.0;
case kSPINFLIP: return 1.0;
case kLOSS: return 1.0; // ??
case kLOSSCS: return CLHEP::barn/units::barn; // ??
case kABSCS: return CLHEP::barn/units::barn; // ??
case kSCATCS: return CLHEP::barn/units::barn; // ??
case kMR_NBTHETA: return 1.0;
case kMR_NBE: return 1.0;
case kMR_RRMS: return 1.0; // ??
case kMR_CORRLEN: return CLHEP::m/units::m; // Length
case kMR_THETAMIN: return 1.0;
case kMR_THETAMAX: return 1.0;
case kMR_EMIN: return CLHEP::keV/units::keV; // Energy
case kMR_EMAX: return CLHEP::keV/units::keV; // Energy
case kMR_ANGNOTHETA: return 1.0;
case kMR_ANGNOPHI: return 1.0;
case kMR_ANGCUT: return 1.0;
default:
break;
}
printout(FATAL,"Geant4Converter", "+++ Cannot convert CONST material property with index: %d", index);
#else
printout(FATAL,"Geant4Converter", "+++ Cannot convert material property with index: %d [Need Geant4 > 10.03]", index);
#endif
Geant4Converter::Geant4Converter(const Detector& description_ref)
: Geant4Mapping(description_ref), checkOverlaps(true) {
this->Geant4Mapping::init();
Markus Frank
committed
m_propagateRegions = true;
Markus Frank
committed
outputLevel = PrintLevel(printLevel() - 1);
Geant4Converter::Geant4Converter(const Detector& description_ref, PrintLevel level)
: Geant4Mapping(description_ref), checkOverlaps(true) {
Markus Frank
committed
m_propagateRegions = true;
Markus Frank
committed
outputLevel = level;
Markus Frank
committed
/// Standard destructor
Geant4Converter::~Geant4Converter() {
Markus Frank
committed
/// Dump element in GDML format to output stream
void* Geant4Converter::handleElement(const string& name, const Atom element) const {
G4Element* g4e = data().g4Elements[element];
if (!g4e) {
g4e = G4Element::GetElement(name, false);
if (!g4e) {
Markus Frank
committed
PrintLevel lvl = debugElements ? ALWAYS : outputLevel;
Markus Frank
committed
double a_conv = (CLHEP::g / CLHEP::mole);
if (element->GetNisotopes() > 0) {
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) {
Markus Frank
committed
g4iso = new G4Isotope(iso->GetName(), iso->GetZ(), iso->GetN(), iso->GetA()*a_conv);
printout(lvl, "Geant4Converter", "++ Created G4 Isotope %s from data: Z=%d N=%d A=%.3f [g/mole]",
iso->GetName(), iso->GetZ(), iso->GetN(), iso->GetA());
}
else {
printout(lvl, "Geant4Converter", "++ Re-used G4 Isotope %s from data: Z=%d N=%d A=%.3f [g/mole]",
Markus Frank
committed
iso->GetName(), iso->GetZ(), iso->GetN(), iso->GetA());
}
g4e->AddIsotope(g4iso, element->GetRelativeAbundance(i));
}
// This adds in Geant4 the natural isotopes, which we normally do not want. We want to steer it outselves.
Markus Frank
committed
g4e = new G4Element(element->GetTitle(), name, element->Z(), element->A()*a_conv);
printout(lvl, "Geant4Converter", "++ Created G4 Isotope %s from data: Z=%d N=%d A=%.3f [g/mole]",
element->GetName(), element->Z(), element->N(), element->A());
#if 0 // Disabled for now!
g4e = new G4Element(name, element->GetTitle(), 1);
G4Isotope* g4iso = G4Isotope::GetIsotope(name, false);
if (!g4iso) {
g4iso = new G4Isotope(name, element->Z(), element->N(), element->A()*a_conv);
printout(lvl, "Geant4Converter", "++ Created G4 Isotope %s from data: Z=%d N=%d A=%.3f [g/mole]",
element->GetName(), element->Z(), element->N(), element->A());
}
else {
printout(lvl, "Geant4Converter", "++ Re-used G4 Isotope %s from data: Z=%d N=%d A=%.3f [g/mole]",
element->GetName(), element->Z(), element->N(), element->A());
}
g4e->AddIsotope(g4iso, 1.0);
#endif
Markus Frank
committed
stringstream str;
str << (*g4e);
printout(lvl, "Geant4Converter", "++ Created G4 %s No.Isotopes:%d",
Markus Frank
committed
str.str().c_str(),element->GetNisotopes());
}
data().g4Elements[element] = g4e;
}
return g4e;
}
/// Dump material in GDML format to output stream
void* Geant4Converter::handleMaterial(const string& name, Material medium) const {
Geant4GeometryInfo& info = data();
G4Material* mat = info.g4Materials[medium];
Markus Frank
committed
PrintLevel lvl = debugMaterials ? ALWAYS : outputLevel;
mat = G4Material::GetMaterial(name, false);
if (!mat) {
Markus Frank
committed
TGeoMaterial* material = medium->GetMaterial();
Markus Frank
committed
G4State state = kStateUndefined;
double density = material->GetDensity() * (CLHEP::gram / CLHEP::cm3);
if (density < 1e-25)
density = 1e-25;
Markus Frank
committed
switch (material->GetState()) {
state = kStateSolid;
break;
state = kStateLiquid;
break;
default:
case TGeoMaterial::kMatStateUndefined:
state = kStateUndefined;
break;
Markus Frank
committed
if (material->IsMixture()) {
double W_total = 0.0;
Markus Frank
committed
TGeoMixture* mix = (TGeoMixture*) material;
int nElements = mix->GetNelements();
Markus Frank
committed
mat = new G4Material(name, density, nElements, state,
material->GetTemperature(), material->GetPressure());
for (int i = 0; i < nElements; ++i) {
A_total += (mix->GetAmixt())[i];
W_total += (mix->GetWmixt())[i];
Markus Frank
committed
}
for (int i = 0; i < nElements; ++i) {
TGeoElement* e = mix->GetElement(i);
G4Element* g4e = (G4Element*) handleElement(e->GetName(), Atom(e));
printout(ERROR, "Material",
"Missing component %s for material %s. A=%f W=%f",
e->GetName(), mix->GetName(), A_total, W_total);
//mat->AddElement(g4e, (mix->GetAmixt())[i] / A_total);
mat->AddElement(g4e, (mix->GetWmixt())[i] / W_total);
Markus Frank
committed
double z = material->GetZ(), a = material->GetA();
Markus Frank
committed
if ( z < 1.0000001 ) z = 1.0;
if ( a < 0.5000001 ) a = 1.0;
Markus Frank
committed
mat = new G4Material(name, z, a, density, state,
material->GetTemperature(), material->GetPressure());
#if ROOT_VERSION_CODE >= ROOT_VERSION(6,17,0)
/// Attach the material properties if any
G4MaterialPropertiesTable* tab = 0;
TListIter propIt(&material->GetProperties());
Markus Frank
committed
for(TObject* obj=propIt.Next(); obj; obj = propIt.Next()) {
TNamed* named = (TNamed*)obj;
TGDMLMatrix* matrix = info.manager->GetGDMLMatrix(named->GetTitle());
Geant4GeometryInfo::PropertyVector* v =
(Geant4GeometryInfo::PropertyVector*)handleMaterialProperties(matrix);
if ( 0 == v ) {
except("Geant4Converter", "++ FAILED to create G4 material %s [Cannot convert property:%s]",
material->GetName(), named->GetName());
Markus Frank
committed
if ( 0 == tab ) {
tab = new G4MaterialPropertiesTable();
mat->SetMaterialPropertiesTable(tab);
}
int idx = tab->GetPropertyIndex(named->GetName(), false);
printout(ERROR, "Geant4Converter", "++ UNKNOWN Geant4 CONST Property: %-20s [IGNORED]", named->GetName());
continue;
}
// We need to convert the property from TGeo units to Geant4 units
auto conv = g4PropertyConversion(idx);
vector<double> bins(v->bins), vals(v->values);
for(size_t i=0, count=bins.size(); i<count; ++i)
bins[i] *= conv.first, vals[i] *= conv.second;
G4MaterialPropertyVector* vec = new G4MaterialPropertyVector(&bins[0], &vals[0], bins.size());
tab->AddProperty(named->GetName(), vec);
Markus Frank
committed
printout(lvl, "Geant4Converter", "++ Property: %-20s [%ld x %ld] -> %s ",
named->GetName(), matrix->GetRows(), matrix->GetCols(), named->GetTitle());
for(size_t i=0, count=v->bins.size(); i<count; ++i)
printout(lvl,named->GetName()," Geant4: %8.3g [MeV] TGeo: %8.3g [GeV] Conversion: %8.3g",
bins[i], v->bins[i], conv.first);
/// Attach the material properties if any
TListIter cpropIt(&material->GetConstProperties());
for(TObject* obj=cpropIt.Next(); obj; obj = cpropIt.Next()) {
Bool_t err = kFALSE;
TNamed* named = (TNamed*)obj;
double v = info.manager->GetProperty(named->GetTitle(),&err);
if ( err != kFALSE ) {
except("Geant4Converter",
"++ FAILED to create G4 material %s "
"[Cannot convert const property: %s]",
material->GetName(), named->GetName());
}
if ( 0 == tab ) {
tab = new G4MaterialPropertiesTable();
mat->SetMaterialPropertiesTable(tab);
}
int idx = tab->GetConstPropertyIndex(named->GetName(), false);
printout(ERROR, "Geant4Converter", "++ UNKNOWN Geant4 CONST Property: %-20s [IGNORED]", named->GetName());
continue;
}
// We need to convert the property from TGeo units to Geant4 units
double conv = g4ConstPropertyConversion(idx);
printout(lvl, "Geant4Converter", "++ CONST Property: %-20s %g ", named->GetName(), v);
tab->AddConstProperty(named->GetName(), v * conv);
auto* ionization = mat->GetIonisation();
stringstream str;
str << (*mat) << endl;
if ( ionization )
str << " log(MEE): " << ionization->GetLogMeanExcEnergy();
else
str << " log(MEE): UNKNOWN";
printout(lvl, "Geant4Converter", "++ Created G4 %s", str.str().c_str());
info.g4Materials[medium] = mat;
void* Geant4Converter::handleSolid(const string& name, const TGeoShape* shape) const {
Markus Frank
committed
G4VSolid* solid = 0;
Markus Frank
committed
PrintLevel lvl = debugShapes ? ALWAYS : outputLevel;
if (0 != (solid = data().g4Solids[shape])) {
Markus Frank
committed
return solid;
}
else if (shape->IsA() == TGeoShapeAssembly::Class()) {
// Assemblies have no corresponding 'shape' in Geant4. Ignore the shape translation.
// It does not harm, since this 'shape' is never accessed afterwards.
data().g4Solids[shape] = solid;
return solid;
Markus Frank
committed
}
else if (shape->IsA() == TGeoBBox::Class()) {
const TGeoBBox* sh = (const TGeoBBox*) shape;
solid = new G4Box(name, sh->GetDX() * CM_2_MM, sh->GetDY() * CM_2_MM, sh->GetDZ() * CM_2_MM);
else if (shape->IsA() == TGeoTube::Class()) {
const TGeoTube* sh = (const TGeoTube*) shape;
solid = new G4Tubs(name, sh->GetRmin() * CM_2_MM, sh->GetRmax() * CM_2_MM, sh->GetDz() * CM_2_MM, 0, 2. * M_PI);
else if (shape->IsA() == TGeoTubeSeg::Class()) {
const TGeoTubeSeg* sh = (const TGeoTubeSeg*) shape;
solid = new G4Tubs(name, sh->GetRmin() * CM_2_MM, sh->GetRmax() * CM_2_MM, sh->GetDz() * CM_2_MM,
sh->GetPhi1() * DEGREE_2_RAD, (sh->GetPhi2()-sh->GetPhi1()) * DEGREE_2_RAD);
else if (shape->IsA() == TGeoEltu::Class()) {
const TGeoEltu* sh = (const TGeoEltu*) shape;
solid = new G4EllipticalTube(name,sh->GetA() * CM_2_MM, sh->GetB() * CM_2_MM, sh->GetDz() * CM_2_MM);
else if (shape->IsA() == TGeoTrd1::Class()) {
const TGeoTrd1* sh = (const TGeoTrd1*) shape;
solid = new G4Trd(name, sh->GetDx1() * CM_2_MM, sh->GetDx2() * CM_2_MM, sh->GetDy() * CM_2_MM, sh->GetDy() * CM_2_MM,
sh->GetDz() * CM_2_MM);
else if (shape->IsA() == TGeoTrd2::Class()) {
const TGeoTrd2* sh = (const TGeoTrd2*) shape;
solid = new G4Trd(name, sh->GetDx1() * CM_2_MM, sh->GetDx2() * CM_2_MM, sh->GetDy1() * CM_2_MM, sh->GetDy2() * CM_2_MM,
sh->GetDz() * CM_2_MM);
else if (shape->IsA() == TGeoHype::Class()) {
const TGeoHype* sh = (const TGeoHype*) shape;
solid = new G4Hype(name, sh->GetRmin() * CM_2_MM, sh->GetRmax() * CM_2_MM,
sh->GetStIn() * DEGREE_2_RAD, sh->GetStOut() * DEGREE_2_RAD,
sh->GetDz() * CM_2_MM);
else if (shape->IsA() == TGeoXtru::Class()) {
const TGeoXtru* sh = (const TGeoXtru*) shape;
size_t nz = sh->GetNz();
vector<G4ExtrudedSolid::ZSection> z;
vector<G4TwoVector> polygon;
z.reserve(nz);
polygon.reserve(nz);
for(size_t i=0; i<nz; ++i) {
z.emplace_back(G4ExtrudedSolid::ZSection(sh->GetZ(i) * CM_2_MM,
{sh->GetXOffset(i), sh->GetYOffset(i)},
sh->GetScale(i)));
polygon.emplace_back(G4TwoVector(sh->GetX(i) * CM_2_MM,sh->GetY(i) * CM_2_MM));
}
solid = new G4ExtrudedSolid(name, polygon, z);
}
else if (shape->IsA() == TGeoPgon::Class()) {
const TGeoPgon* sh = (const TGeoPgon*) shape;
double phi_start = sh->GetPhi1() * DEGREE_2_RAD;
double phi_total = (sh->GetDphi() + sh->GetPhi1()) * DEGREE_2_RAD;
for (Int_t i = 0; i < sh->GetNz(); ++i) {
rmin.emplace_back(sh->GetRmin(i) * CM_2_MM);
rmax.emplace_back(sh->GetRmax(i) * CM_2_MM);
z.emplace_back(sh->GetZ(i) * CM_2_MM);
solid = new G4Polyhedra(name, phi_start, phi_total, sh->GetNedges(), sh->GetNz(), &z[0], &rmin[0], &rmax[0]);
else if (shape->IsA() == TGeoPcon::Class()) {
const TGeoPcon* sh = (const TGeoPcon*) shape;
double phi_start = sh->GetPhi1() * DEGREE_2_RAD;
double phi_total = (sh->GetDphi() + sh->GetPhi1()) * DEGREE_2_RAD;
for (Int_t i = 0; i < sh->GetNz(); ++i) {
rmin.emplace_back(sh->GetRmin(i) * CM_2_MM);
rmax.emplace_back(sh->GetRmax(i) * CM_2_MM);
z.emplace_back(sh->GetZ(i) * CM_2_MM);
solid = new G4Polycone(name, phi_start, phi_total, sh->GetNz(), &z[0], &rmin[0], &rmax[0]);
const TGeoCone* sh = (const TGeoCone*) shape;
solid = new G4Cons(name, sh->GetRmin1() * CM_2_MM, sh->GetRmax1() * CM_2_MM, sh->GetRmin2() * CM_2_MM,
sh->GetRmax2() * CM_2_MM, sh->GetDz() * CM_2_MM, 0.0, 2.*M_PI);
else if (shape->IsA() == TGeoConeSeg::Class()) {
const TGeoConeSeg* sh = (const TGeoConeSeg*) shape;
solid = new G4Cons(name, sh->GetRmin1() * CM_2_MM, sh->GetRmax1() * CM_2_MM,
sh->GetRmin2() * CM_2_MM, sh->GetRmax2() * CM_2_MM,
sh->GetDz() * CM_2_MM,
sh->GetPhi1() * DEGREE_2_RAD, (sh->GetPhi2()-sh->GetPhi1()) * DEGREE_2_RAD);
}
else if (shape->IsA() == TGeoParaboloid::Class()) {
const TGeoParaboloid* sh = (const TGeoParaboloid*) shape;
solid = new G4Paraboloid(name, sh->GetDz() * CM_2_MM, sh->GetRlo() * CM_2_MM, sh->GetRhi() * CM_2_MM);
#if 0 /* Not existent */
else if (shape->IsA() == TGeoEllisoid::Class()) {
const TGeoParaboloid* sh = (const TGeoParaboloid*) shape;
solid = new G4Paraboloid(name, sh->GetDz() * CM_2_MM, sh->GetRlo() * CM_2_MM, sh->GetRhi() * CM_2_MM);
}
#endif
else if (shape->IsA() == TGeoSphere::Class()) {
const TGeoSphere* sh = (const TGeoSphere*) shape;
solid = new G4Sphere(name, sh->GetRmin() * CM_2_MM, sh->GetRmax() * CM_2_MM, sh->GetPhi1() * DEGREE_2_RAD,
sh->GetPhi2() * DEGREE_2_RAD, sh->GetTheta1() * DEGREE_2_RAD, sh->GetTheta2() * DEGREE_2_RAD);
}
else if (shape->IsA() == TGeoTorus::Class()) {
const TGeoTorus* sh = (const TGeoTorus*) shape;
solid = new G4Torus(name, sh->GetRmin() * CM_2_MM, sh->GetRmax() * CM_2_MM, sh->GetR() * CM_2_MM,
sh->GetPhi1() * DEGREE_2_RAD, sh->GetDphi() * DEGREE_2_RAD);
else if (shape->IsA() == TGeoTrap::Class()) {
const TGeoTrap* sh = (const TGeoTrap*) shape;
solid = new G4Trap(name, sh->GetDz() * CM_2_MM, sh->GetTheta() * DEGREE_2_RAD, sh->GetPhi() * DEGREE_2_RAD,
sh->GetH1() * CM_2_MM, sh->GetBl1() * CM_2_MM, sh->GetTl1() * CM_2_MM, sh->GetAlpha1() * DEGREE_2_RAD,
sh->GetH2() * CM_2_MM, sh->GetBl2() * CM_2_MM, sh->GetTl2() * CM_2_MM, sh->GetAlpha2() * DEGREE_2_RAD);
else if (shape->IsA() == TGeoArb8::Class()) {
vector<G4TwoVector> vertices;
TGeoTrap* sh = (TGeoTrap*) shape;
Double_t* vtx_xy = sh->GetVertices();
for ( size_t i=0; i<8; ++i, vtx_xy +=2 )
vertices.emplace_back(G4TwoVector(vtx_xy[0] * CM_2_MM,vtx_xy[1] * CM_2_MM));
solid = new G4GenericTrap(name, sh->GetDz() * CM_2_MM, vertices);
}
else if (shape->IsA() == TGeoScaledShape::Class()) {
TGeoScaledShape* sh = (TGeoScaledShape*) shape;
const double* vals = sh->GetScale()->GetScale();
Solid s_sh(sh->GetShape());
G4VSolid* scaled = (G4VSolid*)handleSolid(s_sh.name(), s_sh.ptr());
solid = new G4ReflectedSolid(scaled->GetName() + "_refl",
scaled, G4Scale3D(vals[0],vals[1],vals[2]));
printout(INFO,"G4Shapes","Converting scaled shape from reflection: %s",sh->GetName());
}
else if (shape->IsA() == TGeoCompositeShape::Class()) {
const TGeoCompositeShape* sh = (const TGeoCompositeShape*) shape;
const TGeoBoolNode* boolean = sh->GetBoolNode();
TGeoBoolNode::EGeoBoolType oper = boolean->GetBooleanOperator();
TGeoMatrix* matrix = boolean->GetRightMatrix();
G4VSolid* left = (G4VSolid*) handleSolid(name + "_left", boolean->GetLeftShape());
G4VSolid* right = (G4VSolid*) handleSolid(name + "_right", boolean->GetRightShape());
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);
//specific case!
//Ellipsoid tag preparing
//if left == TGeoScaledShape AND right == TGeoBBox
// AND if TGeoScaledShape->GetShape == TGeoSphere
TGeoShape* ls = boolean->GetLeftShape();
TGeoShape* rs = boolean->GetRightShape();
if (strcmp(ls->ClassName(), "TGeoScaledShape") == 0 &&
Markus Frank
committed
strcmp(rs->ClassName(), "TGeoBBox") == 0) {
if (strcmp(((TGeoScaledShape *)ls)->GetShape()->ClassName(), "TGeoSphere") == 0) {
if (oper == TGeoBoolNode::kGeoIntersection) {
TGeoScaledShape* lls = (TGeoScaledShape *)ls;
TGeoBBox* rrs = (TGeoBBox*)rs;
Markus Frank
committed
double sx = lls->GetScale()->GetScale()[0];
double sy = lls->GetScale()->GetScale()[1];
Markus Frank
committed
double radius = ((TGeoSphere *)lls->GetShape())->GetRmax();
Markus Frank
committed
double dz = rrs->GetDZ();
double zorig = rrs->GetOrigin()[2];
double zcut2 = dz + zorig;
double zcut1 = 2 * zorig - zcut2;
Markus Frank
committed
solid = new G4Ellipsoid(name,
sx * radius * CM_2_MM,
sy * radius * CM_2_MM,
radius * CM_2_MM,
zcut1 * CM_2_MM,
zcut2 * CM_2_MM);
data().g4Solids[shape] = solid;
return solid;
}
}
if (matrix->IsRotation()) {
MyTransform3D transform(matrix->GetTranslation(),matrix->GetRotationMatrix());
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);
const Double_t *t = matrix->GetTranslation();
G4ThreeVector transform(t[0] * CM_2_MM, t[1] * CM_2_MM, t[2] * CM_2_MM);
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());
Markus Frank
committed
else {
printout(lvl,"Geant4Converter","++ Successessfully converted shape [%p] of type:%s to %s.",
solid,shape->IsA()->GetName(),typeName(typeid(*solid)).c_str());
}
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 {
Geant4GeometryInfo& info = data();
PrintLevel lvl = debugVolumes ? ALWAYS : outputLevel;
Markus Frank
committed
Geant4GeometryMaps::VolumeMap::const_iterator volIt = info.g4Volumes.find(volume);
Volume _v(volume);
if ( _v.testFlagBit(Volume::VETO_SIMU) ) {
printout(lvl, "Geant4Converter", "++ Volume %s not converted [Veto'ed for simulation]",volume->GetName());
return 0;
}
else if (volIt == info.g4Volumes.end() ) {
string n = v->GetName();
TGeoShape* sh = v->GetShape();
G4VSolid* solid = (G4VSolid*) handleSolid(sh->GetName(), sh);
bool assembly = sh->IsA() == TGeoShapeAssembly::Class() || v->IsA() == TGeoVolumeAssembly::Class();
LimitSet lim = _v.limitSet();
Markus Frank
committed
G4UserLimits* user_limits = 0;
if (!user_limits) {
throw runtime_error("G4Cnv::volume[" + name + "]: + FATAL Failed to "
Markus Frank
committed
"access Geant4 user limits.");
VisAttr vis = _v.visAttributes();
G4VisAttributes* vis_attr = 0;
vis_attr = (G4VisAttributes*)handleVis(vis.name(), vis);
Markus Frank
committed
}
Markus Frank
committed
G4Region* region = 0;
if (!region) {
throw runtime_error("G4Cnv::volume[" + name + "]: + FATAL Failed to "
Markus Frank
committed
"access Geant4 region.");
Markus Frank
committed
}
Markus Frank
committed
printout(lvl, "Geant4Converter", "++ Convert Volume %-32s: %p %s/%s assembly:%s",
n.c_str(), v, sh->IsA()->GetName(), v->IsA()->GetName(), yes_no(assembly));
//info.g4AssemblyVolumes[v] = new Geant4AssemblyVolume();
Markus Frank
committed
}
medium = (G4Material*) handleMaterial(med->GetName(), Material(med));
if (!solid) {
throw runtime_error("G4Converter: No Geant4 Solid present for volume:" + n);
Markus Frank
committed
}
if (!medium) {
throw runtime_error("G4Converter: No Geant4 material present for volume:" + n);
Markus Frank
committed
}
printout(lvl, "Geant4Converter", "++ Volume + Apply LIMITS settings:%-24s to volume %s.",
Markus Frank
committed
lim.name(), _v.name());
Markus Frank
committed
}
Markus Frank
committed
G4LogicalVolume* vol = new G4LogicalVolume(solid, medium, n, 0, 0, user_limits);
printout(lvl, "Geant4Converter", "++ Volume + Apply REGION settings: %s to volume %s.",
Markus Frank
committed
reg.name(), _v.name());
Markus Frank
committed
vol->SetRegion(region);
region->AddRootLogicalVolume(vol);
Markus Frank
committed
vol->SetVisAttributes(vis_attr);
}
info.g4Volumes[v] = vol;
printout(lvl, "Geant4Converter", "++ Volume + %s converted: %p ---> G4: %p", n.c_str(), v, vol);
return nullptr;
/// Dump logical volume in GDML format to output stream
void* Geant4Converter::collectVolume(const string& /* name */, const TGeoVolume* volume) const {
Geant4GeometryInfo& info = data();
Region reg = _v.region();
LimitSet lim = _v.limitSet();
SensitiveDetector det = _v.sensitiveDetector();
/// Dump volume placement in GDML format to output stream
void* Geant4Converter::handleAssembly(const string& name, const TGeoNode* node) const {
TGeoVolume* mot_vol = node->GetVolume();
PrintLevel lvl = debugVolumes ? ALWAYS : outputLevel;
if ( mot_vol->IsA() != TGeoVolumeAssembly::Class() ) {
Volume _v(mot_vol);
if ( _v.testFlagBit(Volume::VETO_SIMU) ) {
printout(lvl, "Geant4Converter", "++ AssemblyNode %s not converted [Veto'ed for simulation]",node->GetName());
return 0;
}
Geant4GeometryInfo& info = data();
Geant4AssemblyVolume* g4 = info.g4AssemblyVolumes[node];
if ( !g4 ) {
g4 = new Geant4AssemblyVolume();
for(Int_t i=0; i < mot_vol->GetNdaughters(); ++i) {
TGeoNode* d = mot_vol->GetNode(i);
TGeoVolume* dau_vol = d->GetVolume();
TGeoMatrix* tr = d->GetMatrix();
MyTransform3D transform(tr->GetTranslation(),tr->IsRotation() ? tr->GetRotationMatrix() : s_identity_rot);
if ( dau_vol->IsA() == TGeoVolumeAssembly::Class() ) {
Markus Frank
committed
Geant4GeometryMaps::AssemblyMap::iterator assIt = info.g4AssemblyVolumes.find(d);
Markus Frank
committed
if ( assIt == info.g4AssemblyVolumes.end() ) {
printout(FATAL, "Geant4Converter", "+++ Invalid child assembly at %s : %d parent: %s child:%s",
__FILE__, __LINE__, name.c_str(), d->GetName());
Markus Frank
committed
}
g4->placeAssembly(d,(*assIt).second,transform);
printout(lvl, "Geant4Converter", "+++ Assembly: AddPlacedAssembly : dau:%s "
Markus Frank
committed
"to mother %s Tr:x=%8.3f y=%8.3f z=%8.3f",
dau_vol->GetName(), mot_vol->GetName(),
transform.dx(), transform.dy(), transform.dz());
Markus Frank
committed
Geant4GeometryMaps::VolumeMap::iterator volIt = info.g4Volumes.find(dau_vol);
Markus Frank
committed
if ( volIt == info.g4Volumes.end() ) {
printout(FATAL,"Geant4Converter", "+++ Invalid child volume at %s : %d parent: %s child:%s",
Markus Frank
committed
__FILE__, __LINE__, name.c_str(), d->GetName());
except("Geant4Converter", "+++ Invalid child volume at %s : %d parent: %s child:%s",
__FILE__, __LINE__, name.c_str(), d->GetName());
Markus Frank
committed
}
g4->placeVolume(d,(*volIt).second, transform);
printout(lvl, "Geant4Converter", "+++ Assembly: AddPlacedVolume : dau:%s "
Markus Frank
committed
"to mother %s Tr:x=%8.3f y=%8.3f z=%8.3f",
dau_vol->GetName(), mot_vol->GetName(),
transform.dx(), transform.dy(), transform.dz());
info.g4AssemblyVolumes[node] = g4;
/// Dump volume placement in GDML format to output stream
void* Geant4Converter::handlePlacement(const string& name, const TGeoNode* node) const {
Geant4GeometryInfo& info = data();
PrintLevel lvl = debugPlacements ? ALWAYS : outputLevel;
Markus Frank
committed
Geant4GeometryMaps::PlacementMap::const_iterator g4it = info.g4Placements.find(node);
G4VPhysicalVolume* g4 = (g4it == info.g4Placements.end()) ? 0 : (*g4it).second;
TGeoVolume* vol = node->GetVolume();
Volume _v(vol);
if ( _v.testFlagBit(Volume::VETO_SIMU) ) {
printout(lvl, "Geant4Converter", "++ Placement %s not converted [Veto'ed for simulation]",node->GetName());
return 0;
}
if (!g4) {
TGeoVolume* mot_vol = node->GetMotherVolume();
TGeoMatrix* tr = node->GetMatrix();
if (!tr) {
except("Geant4Converter", "++ Attempt to handle placement without transformation:%p %s of type %s vol:%p", node,
node->GetName(), node->IsA()->GetName(), vol);
except("Geant4Converter", "++ Unknown G4 volume:%p %s of type %s ptr:%p", node, node->GetName(),
node->IsA()->GetName(), vol);
}
else {
int copy = node->GetNumber();
bool node_is_assembly = vol->IsA() == TGeoVolumeAssembly::Class();
bool mother_is_assembly = mot_vol ? mot_vol->IsA() == TGeoVolumeAssembly::Class() : false;
MyTransform3D transform(tr->GetTranslation(),tr->IsRotation() ? tr->GetRotationMatrix() : s_identity_rot);
Markus Frank
committed
Geant4GeometryMaps::VolumeMap::const_iterator volIt = info.g4Volumes.find(mot_vol);
if ( mother_is_assembly ) {
Markus Frank
committed
//
// Mother is an assembly:
// Nothing to do here, because:
// -- placed volumes were already added before in "handleAssembly"
// -- imprint cannot be made, because this requires a logical volume as a mother
//
printout(lvl, "Geant4Converter", "+++ Assembly: **** : dau:%s "
Markus Frank
committed
"to mother %s Tr:x=%8.3f y=%8.3f z=%8.3f",
vol->GetName(), mot_vol->GetName(),
transform.dx(), transform.dy(), transform.dz());
Markus Frank
committed
}
else if ( node_is_assembly ) {
Markus Frank
committed
//
// Node is an assembly:
// Imprint the assembly. The mother MUST already be transformed.
//
printout(lvl, "Geant4Converter", "+++ Assembly: makeImprint: dau:%s in mother %s "
Markus Frank
committed
"Tr:x=%8.3f y=%8.3f z=%8.3f",
node->GetName(), mot_vol ? mot_vol->GetName() : "<unknown>",
Markus Frank
committed
transform.dx(), transform.dy(), transform.dz());
Geant4AssemblyVolume* ass = (Geant4AssemblyVolume*)info.g4AssemblyVolumes[node];
Geant4AssemblyVolume::Chain chain;
chain.emplace_back(node);
Markus Frank
committed
ass->imprint(info,node,chain,ass,(*volIt).second, transform, copy, checkOverlaps);
Markus Frank
committed
return 0;
else if ( node != info.manager->GetTopNode() && volIt == info.g4Volumes.end() ) {
Markus Frank
committed
throw logic_error("Geant4Converter: Invalid mother volume found!");