Newer
Older
// $Id:$
//====================================================================
// AIDA Detector description implementation for LCD
//--------------------------------------------------------------------
//
// Author : M.Frank
//
//====================================================================
#include "DD4hep/LCDD.h"
#include "DDG4/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 "G4ProductionCuts.hh"
#include "G4VUserRegionInformation.hh"
// Geant4 include files
#include "G4Element.hh"
#include "G4Box.hh"
#include "G4Tubs.hh"
#include "G4Trd.hh"
#include "G4Paraboloid.hh"
#include "G4Polycone.hh"
#include "G4Polyhedra.hh"
#include "G4Sphere.hh"
#include "G4Torus.hh"
#include "G4UnionSolid.hh"
#include "G4SubtractionSolid.hh"
#include "G4IntersectionSolid.hh"
#include "G4Region.hh"
#include "G4UserLimits.hh"
#include "G4VSensitiveDetector.hh"
#include "G4LogicalVolume.hh"
#include "G4Material.hh"
#include "G4Element.hh"
#include "G4Isotope.hh"
#include "G4Transform3D.hh"
#include "G4ThreeVector.hh"
#include "G4PVPlacement.hh"
using namespace DD4hep::Simulation;
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); }
};
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;
}
class G4UserRegionInformation : public G4VUserRegionInformation {
public:
Region region;
double threshold;
bool storeSecondaries;
G4UserRegionInformation() : threshold(0.0), storeSecondaries(false) {}
virtual ~G4UserRegionInformation() {}
virtual void Print() const {
if ( region.isValid() ) cout << "Region:" << region.name() << endl;
}
};
}
/// 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 {
G4double z = element->Z();
G4double aeff = element->A()*(g/mole);
cout << name << " " << g/mole << " " << element->A() << " " << aeff/(g/mole)
<< " " << (aeff/(g/mole) - z) << endl;
g4e = new G4Element(element->GetTitle(),name,z,aeff);
}
}
data().g4Elements[element] = g4e;
}
return g4e;
}
/// Dump material in GDML format to output stream
void* Geant4Converter::handleMaterial(const string& name, const TGeoMedium* medium) const {
G4Material* mat = data().g4Materials[medium];
if ( !mat ) {
mat = G4Material::GetMaterial(name,false);
if ( !mat ) {
int nElements = 0;
const char* opt = "----";
TGeoMaterial* m = medium->GetMaterial();
if ( m->IsMixture() ) {
double A_total = 0.0;
nElements = mix->GetNelements();
mat = new G4Material(name,mix->GetDensity(),nElements);
for(int i=0; i<nElements; ++i)
A_total += (mix->GetAmixt())[i];
for(int i=0; i<nElements; ++i) {
TGeoElement* e = mix->GetElement(i);
G4Element* g4e = (G4Element*)handleElement(e->GetName(),e);
if ( !g4e ) {
cout << "ERROR: Missing component " << e->GetName()
<< " for material " << mix->GetName() << endl;
}
mat->AddElement(g4e,(mix->GetAmixt())[i]/A_total);
opt = "mixture";
}
else {
mat = new G4Material(name,m->GetZ(),m->GetA(),m->GetDensity());
opt = "raw_material";
cout << "Created G4 Material " << m->GetName() << " as " << opt << " " << nElements << endl;
}
data().g4Materials[medium] = mat;
}
return mat;
}
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() ) {
const TGeoPgon* s = (const TGeoPgon*)shape;
double phi_start = s->GetPhi1()*DEGREE_2_RAD;
double phi_total = (s->GetDphi()+s->GetPhi1())*DEGREE_2_RAD;
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 G4Polyhedra(name,phi_start,phi_total,s->GetNedges()-1,s->GetNz(),&z[0],&rmin[0],&rmax[0]);
}
else if ( shape->IsA() == TGeoPcon::Class() ) {
const TGeoPcon* s = (const TGeoPcon*)shape;
double phi_start = s->GetPhi1()*DEGREE_2_RAD;
double phi_total = (s->GetDphi()+s->GetPhi1())*DEGREE_2_RAD;
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,phi_start,phi_total,s->GetNz(),&z[0],&rmin[0],&rmax[0]);
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
}
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 {
G4LogicalVolume* vol = data().g4Volumes[volume];
const TGeoVolume* v = volume;
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;
// Here we collect all information we need later to handle to volume internals
Volume _v = Ref_t(v);
Region reg = _v.region();
LimitSet lim = _v.limitSet();
SensitiveDetector det = _v.sensitiveDetector();
if ( lim.isValid() ) data().limits[lim.ptr()].insert(v);
if ( reg.isValid() ) data().regions[reg.ptr()].insert(v);
if ( det.isValid() ) data().sensitives[det.ptr()].insert(v);
//cout << "Converted logical volume [" << n << "]:" << v.ptr() << " ---> G4 " << vol << endl;
}
return vol;
}
/// Dump volume placement in GDML format to output stream
void* Geant4Converter::handlePlacement(const string& name, const TGeoNode* node) const {
G4PVPlacement* g4 = data().g4Placements[node];
if ( !g4 ) {
TGeoMatrix* trafo = node->GetMatrix();
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()];
const Value<TGeoNodeMatrix,PlacedVolume::Object>* obj =
dynamic_cast<const Value<TGeoNodeMatrix,PlacedVolume::Object>* >(node);
if ( 0 == g4vol ) {
cout << "FATAL: Unknown G4 volume:" << (void*)obj << " " << obj->GetName() << endl;
}
else 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]);
g4 = 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]);
g4 = 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
// cout << "Created volume placement:" << name << " No:" << copy << endl;
}
}
else {
cout << "Attempt to DOUBLE-place physical volume:" << name << " No:" << node->GetNumber() << endl;
}
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
return g4;
}
/// Convert the geometry type region into the corresponding Geant4 object(s).
void* Geant4Converter::handleRegion(const TNamed* region, const set<const TGeoVolume*>& volumes) const {
G4Region* g4 = data().g4Regions[region];
if ( !g4 ) {
Region r = Ref_t(region);
g4 = new G4Region(region->GetName());
// set production cut
G4ProductionCuts* cuts = new G4ProductionCuts();
cuts->SetProductionCut(r.cut());
g4->SetProductionCuts(cuts);
// create region info with storeSecondaries flag
G4UserRegionInformation* info = new G4UserRegionInformation();
info->region = r;
info->threshold = r.threshold();
info->storeSecondaries = r.storeSecondaries();
g4->SetUserInformation(info);
vector<string>& limits = r.limits();
for(vector<string>::const_iterator i=limits.begin(); i!=limits.end(); ++i) {
const string& nam = *i;
LimitSet ls = m_lcdd.limitSet(nam);
if ( ls.isValid() ) {
bool found = false;
const LimitMap& lm = data().g4Limits;
for(LimitMap::const_iterator j=lm.begin(); j!=lm.end();++j) {
if ( nam == (*j).first->GetName() ) {
g4->SetUserLimits((*j).second);
found = true;
break;
}
}
if ( found ) continue;
}
throw runtime_error("G4Region: Failed to resolve user limitset:"+*i);
}
data().g4Regions[region] = g4;
}
return g4;
}
/// Convert the geometry type LimitSet into the corresponding Geant4 object(s).
void* Geant4Converter::handleLimitSet(const TNamed* limitset, const set<const TGeoVolume*>& volumes) const {
G4UserLimits* g4 = data().g4Limits[limitset];
if ( !g4 ) {
LimitSet ls = Ref_t(limitset);
g4 = new G4UserLimits(limitset->GetName());
const LimitSet::Object& obj = ls.limits();
for(LimitSet::Object::const_iterator i=obj.begin(); i!=obj.end(); ++i) {
const Limit& l = *i;
if ( l.name == "step_length_max" )
g4->SetMaxAllowedStep(l.value);
else if ( l.name == "track_length_max" )
g4->SetMaxAllowedStep(l.value);
else if ( l.name == "time_max" )
g4->SetUserMaxTime(l.value);
else if ( l.name == "ekin_min" )
g4->SetUserMinEkine(l.value);
else if ( l.name == "range_min" )
g4->SetUserMinRange(l.value);
else
throw runtime_error("Unknown Geant4 user limit: "+l.toString());
}
data().g4Limits[limitset] = g4;
}
return g4;
}
/// Convert the geometry type SensitiveDetector into the corresponding Geant4 object(s).
void* Geant4Converter::handleSensitive(const TNamed* sens_det, const set<const TGeoVolume*>& volumes) const {
SensitiveDetector sd = Ref_t(sens_det);
/// LCDD: Sensitive detector type == item tag name
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
string type = sd.type(), name = sd.name();
g4 = ROOT::Reflex::PluginService::Create<Geant4SensitiveDetector*>(type,name,&m_lcdd);
if ( !g4 ) {
throw runtime_error("Geant4Converter<SensitiveDetector>: FATAL Failed to "
"create Geant4 sensitive detector "+name+" of type "+type+".");
}
g4->defineCollection(sd.hitsCollection());
ConstVolumeSet& volset = info.sensitives[sens_det];
for(ConstVolumeSet::iterator i=volset.begin(); i!=volset.end();++i) {
std::map<const TGeoVolume*, G4LogicalVolume*>::iterator v = info.g4Volumes.find(*i);
if ( v == info.g4Volumes.end() ) {
throw runtime_error("Geant4Converter<SensitiveDetector>: FATAL Missing converted "
"Geant 4 logical volume");
}
(*v).second->SetSensitiveDetector(g4);
}
if ( sd.verbose() ) {
cout << "Geant4Converter<SensitiveDetector> +" << setw(18) << left << sd.name()
<< setw(20) << left << " ["+sd.type()+"]"
<< " Hits:" << setw(16) << left << sd.hitsCollection() << endl;
cout << " | "
<< "Cutoff:" << setw(6) << left << sd.energyCutoff()
<< setw(5) << right << volset.size() << " volumes ";
if ( sd.region().isValid() ) cout << " Region:" << setw(12) << left << sd.region().name();
if ( sd.limits().isValid() ) cout << " Limits:" << setw(12) << left << sd.limits().name();
cout << "." << endl;
}
data().g4SensDets[sens_det] = g4;
}
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)->GetName(),*i);
template <typename O, typename C, typename F> void handleMap(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);
}
void Geant4Converter::create(DetElement top) {
G4GeometryInfo& geo = *(m_dataPtr=new G4GeometryInfo);
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);
cout << "++ Handled " << geo.materials.size() << " materials." << endl;
handle(this, geo.solids, &Geant4Converter::handleSolid);
cout << "++ Handled " << geo.solids.size() << " solids." << endl;
handle(this, geo.volumes, &Geant4Converter::handleVolume);
cout << "++ Handled " << geo.volumes.size() << " volumes." << endl;
handleMap(this, geo.limits, &Geant4Converter::handleLimitSet);
cout << "++ Handled " << geo.limits.size() << " limit sets." << endl;
handleMap(this, geo.regions, &Geant4Converter::handleRegion);
cout << "++ Handled " << geo.regions.size() << " regions." << endl;
handleMap(this, geo.sensitives, &Geant4Converter::handleSensitive);
cout << "++ Handled " << geo.sensitives.size() << " sensitive detectors." << endl;
// Now place all this stuff appropriately
for(Data::const_reverse_iterator i=m_data->rbegin(); i != m_data->rend(); ++i)
handle(this, (*i).second, &Geant4Converter::handlePlacement);
}