29 #include <TGeoShape.h>
32 #include <TGeoBoolNode.h>
33 #include <TGeoCompositeShape.h>
37 #include <TGeoMatrix.h>
38 #include <TGeoParaboloid.h>
42 #include <TGeoShapeAssembly.h>
43 #include <TGeoSphere.h>
44 #include <TGeoTorus.h>
48 #include <TGeoScaledShape.h>
66 XYZRotation getXYZangles(
const Double_t* r) {
67 Double_t cosb = std::sqrt(r[0]*r[0] + r[1]*r[1]);
69 return XYZRotation(atan2(r[5], r[8]), atan2(-r[2], cosb), atan2(r[1], r[0]));
71 return XYZRotation(atan2(-r[7], r[4]),atan2(-r[2], cosb),0);
75 XYZRotation getXYZangles(
const Double_t* rotationMatrix) {
78 const Double_t *r = rotationMatrix;
79 Double_t cosb = TMath::Sqrt(r[0] * r[0] + r[1] * r[1]);
81 a = TMath::ATan2(r[5], r[8]) * rad;
82 b = TMath::ATan2(-r[2], cosb) * rad;
83 c = TMath::ATan2(r[1], r[0]) * rad;
86 a = TMath::ATan2(-r[7], r[4]) * rad;
87 b = TMath::ATan2(-r[2], cosb) * rad;
90 XYZRotation rr(a, b, c);
91 std::cout <<
" X:" << a <<
" " << rr.X() <<
" Y:" << b <<
" " << rr.Y() <<
" Z:" << c <<
" " << rr.Z()
92 <<
" lx:" << r[0] <<
" ly:" << r[4] <<
" lz:" << r[8] << std::endl;
93 return XYZRotation(a, b, c);
97 bool is_volume(
const TGeoVolume* volume) {
102 return node.
data() != 0;
105 std::string genName(
const std::string& n) {
return n; }
106 std::string genName(
const std::string& n,
const void* ptr) {
107 std::string nn = genName(n);
109 ::snprintf(text,
sizeof(text),
"%p",ptr);
117 std::map<std::string, const TNamed*>::const_iterator i = _m.find(name);
119 const char* isa = _n ? _n->IsA()->GetName() : (*i).second ? (*i).second->IsA()->GetName() :
"Unknown";
120 std::cout << isa <<
"(position): duplicate entry with name:" << name <<
" " << (
void*) _n <<
" " << (
void*) (*i).second << std::endl;
122 _m.insert(make_pair(name, _n));
141 int Z = element->Z();
142 double A = element->A();
149 e.
setAttr(
_U(formula), element->GetName());
153 atom.
setAttr(
_U(value), A>0.99 ? A : 1.00794 );
165 TGeoMaterial* geo_mat = medium->GetMaterial();
166 double d = geo_mat->GetDensity();
167 if (d < 1e-10) d = 1e-10;
169 mat.
setAttr(
_U(name), medium->GetName());
177 if (geo_mat->IsMixture()) {
178 TGeoMixture *
mix = (TGeoMixture*)geo_mat;
179 const double *wmix =
mix->GetWmixt();
180 const int *nmix =
mix->GetNmixt();
182 for (
int i = 0, n =
mix->GetNelements(); i < n; i++) {
183 TGeoElement *elt =
mix->GetElement(i);
187 for (
int i = 0, n =
mix->GetNelements(); i < n; i++) {
188 TGeoElement *elt =
mix->GetElement(i);
201 else if ( name !=
"dummy" ) {
204 TGeoElement *elt = geo_mat->GetElement(0);
205 printout(INFO,
"++ Converting non mixing material: %s",name.c_str());
212 atom.
setAttr(
_U(value), geo_mat->GetA() );
223 SolidMap::iterator sit = geo.
xmlSolids.find(shape);
231 return (*sit).second;
233 else if (shape->IsA() == TGeoShapeAssembly::Class()) {
240 TClass* isa = shape->IsA();
241 std::string shape_name = shape->GetName();
243 if (isa == TGeoBBox::Class()) {
244 const TGeoBBox* sh = (
const TGeoBBox*) shape;
252 else if (isa == TGeoTube::Class()) {
253 const TGeoTube* sh = (
const TGeoTube*) shape;
264 else if (isa == TGeoTubeSeg::Class()) {
265 const TGeoTubeSeg* sh = (
const TGeoTubeSeg*) shape;
271 solid.
setAttr(
_U(startphi), sh->GetPhi1());
272 solid.
setAttr(
_U(deltaphi), sh->GetPhi2());
276 else if (isa == TGeoEltu::Class()) {
277 const TGeoEltu* sh = (
const TGeoEltu*) shape;
285 else if (isa == TGeoTrd1::Class()) {
286 const TGeoTrd1* sh = (
const TGeoTrd1*) shape;
296 else if (isa == TGeoTrd2::Class()) {
297 const TGeoTrd2* sh = (
const TGeoTrd2*) shape;
307 else if (isa == TGeoHype::Class()) {
308 const TGeoHype* sh = (
const TGeoHype*) shape;
314 solid.
setAttr(
_U(outst), sh->GetStOut());
319 else if (isa == TGeoPgon::Class()) {
320 const TGeoPgon* sh = (
const TGeoPgon*) shape;
323 solid.
setAttr(
_U(startphi), sh->GetPhi1());
324 solid.
setAttr(
_U(deltaphi), sh->GetDphi());
325 solid.
setAttr(
_U(numsides), sh->GetNedges());
328 for (Int_t i = 0; i < sh->GetNz(); ++i) {
331 zplane.
setAttr(
_U(rmin), sh->GetRmin(i));
332 zplane.
setAttr(
_U(rmax), sh->GetRmax(i));
336 else if (isa == TGeoPcon::Class()) {
337 const TGeoPcon* sh = (
const TGeoPcon*) shape;
340 solid.
setAttr(
_U(startphi), sh->GetPhi1());
341 solid.
setAttr(
_U(deltaphi), sh->GetDphi());
344 for (Int_t i = 0; i < sh->GetNz(); ++i) {
347 zplane.
setAttr(
_U(rmin), sh->GetRmin(i));
348 zplane.
setAttr(
_U(rmax), sh->GetRmax(i));
353 else if (isa == TGeoCone::Class()) {
354 const TGeoCone* sh = (
const TGeoCone*) shape;
358 solid.
setAttr(
_U(rmin1), sh->GetRmin1());
359 solid.
setAttr(
_U(rmax1), sh->GetRmax1());
360 solid.
setAttr(
_U(rmin2), sh->GetRmin2());
361 solid.
setAttr(
_U(rmax2), sh->GetRmax2());
367 else if (isa == TGeoConeSeg::Class()) {
372 solid.
setAttr(
_U(rmin1), sh->GetRmin1());
373 solid.
setAttr(
_U(rmin2), sh->GetRmin2());
374 solid.
setAttr(
_U(rmax1), sh->GetRmax1());
375 solid.
setAttr(
_U(rmax2), sh->GetRmax2());
376 solid.
setAttr(
_U(startphi), sh->GetPhi1());
377 solid.
setAttr(
_U(deltaphi), sh->GetPhi2() - sh->GetPhi1());
381 else if (isa == TGeoParaboloid::Class()) {
382 const TGeoParaboloid* sh = (
const TGeoParaboloid*) shape;
391 else if (isa == TGeoEllipsoid::Class()) {
392 const TGeoEllipsoid* sh = (
const TGeoEllipsoid*) shape;
397 else if (isa == TGeoSphere::Class()) {
398 const TGeoSphere* sh = (
const TGeoSphere*) shape;
403 solid.
setAttr(
_U(startphi), sh->GetPhi1());
404 solid.
setAttr(
_U(deltaphi), (sh->GetPhi2() - sh->GetPhi1()));
405 solid.
setAttr(
_U(starttheta), sh->GetTheta1());
406 solid.
setAttr(
_U(deltatheta), (sh->GetTheta2() - sh->GetTheta1()));
410 else if (isa == TGeoTorus::Class()) {
411 const TGeoTorus* sh = (
const TGeoTorus*) shape;
417 solid.
setAttr(
_U(startphi), sh->GetPhi1());
418 solid.
setAttr(
_U(deltaphi), sh->GetDphi());
422 else if (isa == TGeoTrap::Class()) {
423 const TGeoTrap* sh = (
const TGeoTrap*) shape;
433 solid.
setAttr(
_U(alpha1), sh->GetAlpha1());
434 solid.
setAttr(
_U(alpha2), sh->GetAlpha2());
435 solid.
setAttr(
_U(theta), sh->GetTheta());
440 else if (isa == TGeoPara::Class()) {
441 const TGeoPara* sh = (
const TGeoPara*) shape;
447 solid.
setAttr(
_U(alpha), sh->GetAlpha());
448 solid.
setAttr(
_U(theta), sh->GetTheta());
453 else if (isa == TGeoArb8::Class()) {
454 TGeoArb8* sh = (TGeoArb8*) shape;
455 const double* vtx = sh->GetVertices();
477 else if (isa == TGeoScaledShape::Class()) {
478 TGeoScaledShape* sh = (TGeoScaledShape*) shape;
479 const double* vals = sh->GetScale()->GetScale();
480 Solid s_sh(sh->GetShape());
491 else if (isa == TGeoCompositeShape::Class() ||
492 isa == TGeoUnion::Class() ||
493 isa == TGeoIntersection::Class() ||
494 isa == TGeoSubtraction::Class() ) {
495 const TGeoCompositeShape* sh = (
const TGeoCompositeShape*) shape;
496 const TGeoBoolNode*
boolean = sh->GetBoolNode();
497 TGeoBoolNode::EGeoBoolType oper =
boolean->GetBooleanOperator();
498 TGeoMatrix* rm =
boolean->GetRightMatrix();
499 TGeoMatrix* lm =
boolean->GetLeftMatrix();
500 TGeoShape* ls =
boolean->GetLeftShape();
501 TGeoShape* rs =
boolean->GetRightShape();
504 xml_h first_solid(0), second_solid(0);
506 throw std::runtime_error(
"G4Converter: No left Detector Solid present for composite shape:" + name);
509 throw std::runtime_error(
"G4Converter: No right Detector Solid present for composite shape:" + name);
516 if (strcmp(ls->ClassName(),
"TGeoScaledShape") == 0 &&
517 strcmp(rs->ClassName(),
"TGeoBBox") == 0) {
518 if (strcmp(((TGeoScaledShape *)ls)->GetShape()->ClassName(),
"TGeoSphere") == 0) {
519 if (oper == TGeoBoolNode::kGeoIntersection) {
520 TGeoScaledShape* lls = (TGeoScaledShape *)ls;
521 TGeoBBox* rrs = (TGeoBBox*)rs;
524 double sx = lls->GetScale()->GetScale()[0];
525 double sy = lls->GetScale()->GetScale()[1];
526 double radius = ((TGeoSphere *)lls->GetShape())->GetRmax();
527 double dz = rrs->GetDZ();
528 double zorig = rrs->GetOrigin()[2];
529 double zcut2 = dz + zorig;
530 double zcut1 = 2 * zorig - zcut2;
542 if ( oper == TGeoBoolNode::kGeoSubtraction )
544 else if ( oper == TGeoBoolNode::kGeoUnion )
546 else if ( oper == TGeoBoolNode::kGeoIntersection )
550 std::string lnam = left.
attr<std::string>(
_U(name));
551 std::string rnam = right.
attr<std::string>(
_U(name));
556 first_solid.setAttr(
_U(ref), lnam);
557 const double *tr = lm->GetTranslation();
559 if ((tr[0] != 0.0) || (tr[1] != 0.0) || (tr[2] != 0.0)) {
561 obj.
setAttr(
_U(name), name+
"_"+lnam+
"_pos");
567 if (lm->IsRotation()) {
568 TGeoMatrix
const & linv = lm->Inverse();
569 XYZRotation rot = getXYZangles(linv.GetRotationMatrix());
570 if ((rot.X() != 0.0) || (rot.Y() != 0.0) || (rot.Z() != 0.0)) {
572 obj.
setAttr(
_U(name), name+
"_"+lnam+
"_rot");
579 tr = rm->GetTranslation();
582 if ((tr[0] != 0.0) || (tr[1] != 0.0) || (tr[2] != 0.0)) {
586 if (rm->IsRotation()) {
587 TGeoMatrix
const & rinv = rm->Inverse();
588 XYZRotation rot = getXYZangles(rinv.GetRotationMatrix());
589 if ((rot.X() != 0.0) || (rot.Y() != 0.0) || (rot.Z() != 0.0)) {
596 std::string err =
"Failed to handle unknown solid shape:" + name +
" of type " + std::string(shape->IsA()->GetName());
597 throw std::runtime_error(err);
608 const double* tr = trafo->GetTranslation();
609 if (tr[0] != 0.0 || tr[1] != 0.0 || tr[2] != 0.0) {
610 std::string gen_name = genName(name,trafo);
642 XYZRotation r = getXYZangles(trafo->GetRotationMatrix());
643 if (!(r.X() == 0.0 && r.Y() == 0.0 && r.Z() == 0.0)) {
644 std::string gen_name = genName(name,trafo);
676 const TGeoVolume*
v = volume;
678 std::string n = genName(
v->GetName(),
v);
679 TGeoMedium* medium =
v->GetMedium();
680 TGeoShape* sh =
v->GetShape();
684 if (
v->IsAssembly()) {
690 throw std::runtime_error(
"G4Converter: No Detector Solid present for volume:" + n);
692 throw std::runtime_error(
"G4Converter: No Detector material present for volume:" + n);
697 std::string mat_name = medium->GetName();
705 const TObjArray* dau =
const_cast<TGeoVolume*
>(
v)->GetNodes();
706 if (dau && dau->GetEntries() > 0) {
707 for (Int_t i = 0, n_dau = dau->GetEntries(); i < n_dau; ++i) {
708 TGeoNode* node =
reinterpret_cast<TGeoNode*
>(dau->At(i));
743 const TGeoVolume*
v = volume;
745 if (is_volume(volume)) {
762 if ( is_volume(volume) ) {
775 printout(WARNING,
"LCDDConverter",
"++ CollectVolume: Skip volume: %s",volume->GetName());
783 std::stringstream str;
784 str <<
"++ CheckVolumes: Volume " << n <<
" ";
785 if (is_volume(
v.ptr())) {
789 str <<
"of " << sd.
name() <<
" ";
792 str <<
"with VisAttrs " << vis.
name() <<
" ";
795 str <<
"has duplicate entries." << std::endl;
796 printout(ERROR,
"LCDDConverter",str.str().c_str());
807 TGeoMatrix* matrix = node->GetMatrix();
808 TGeoVolume* volume = node->GetVolume();
820 if ( matrix->IsRotation() ) {
826 if (is_placement(node)) {
828 for (PlacedVolume::VolIDs::const_iterator i = ids.begin(); i != ids.end(); ++i) {
830 pvid.
setAttr(
_U(field_name), (*i).first);
839 std::cout <<
"Attempt to DOUBLE-place physical volume:" << name <<
" No:" << node->GetNumber() << std::endl;
867 const std::set<Limit>& obj = lim.
limits();
868 for (std::set<Limit>::const_iterator i = obj.begin(); i != obj.end(); ++i) {
887 std::string typ = seg.
type();
890 for (_P::const_iterator i = p.begin(); i != p.end(); ++i) {
891 const _P::value_type&
v = *i;
892 if (
v->
name() ==
"lunit") {
895 v->value() ==
_toDouble(
"nanometer") ?
"namometer" :
"??";
901 double value =
_toDouble(
v->value()) / dd4hep::mm;
945 int length = 0, start = 0;
948 id.setAttr(
_U(name), name);
950 for (
const auto& i : fm ) {
955 length = f.second<0 ? -f.second : f.second;
956 idfield.
setAttr(
_U(
signed),f.second<0 ?
true :
false);
969 id.setAttr(
_U(length), length + start);
980 float red = 0, green = 0, blue = 0;
995 vis.
setAttr(
_U(drawing_style),
"wireframe");
998 attr.
rgb(red, green, blue);
1015 std::string type = f->GetTitle();
1018 fld = PluginService::Create<NamedObject*>(type +
"_Convert2Detector", &
m_detDesc, &field, &fld);
1019 printout(ALWAYS,
"LCDDConverter",
"++ %s electromagnetic field:%s of type %s",
1020 (fld.
isValid() ?
"Converted" :
"FAILED to convert "), f->GetName(), type.c_str());
1023 PluginService::Create<NamedObject*>(type +
"_Convert2Detector", &
m_detDesc, &field, &fld);
1024 except(
"LCDDConverter",
"Failed to locate plugin to convert electromagnetic field:"
1025 + std::string(f->GetName()) +
" of type " + type +
". "
1035 std::map<std::string, std::string> processors;
1036 static int s_idd = 9999999;
1038 for (Detector::Properties::const_iterator i = prp.begin(); i != prp.end(); ++i) {
1039 const std::string& nam = (*i).first;
1041 if (nam.substr(0, 6) ==
"geant4") {
1042 Detector::PropertyValues::const_iterator id_it = vals.find(
"id");
1043 if (id_it != vals.end()) {
1044 id = (*id_it).second;
1048 ::snprintf(text,
sizeof(text),
"%d", ++s_idd);
1051 processors.insert(make_pair(
id, nam));
1054 for (std::map<std::string, std::string>::const_iterator i = processors.begin(); i != processors.end(); ++i) {
1056 std::string nam = (*i).second;
1058 auto iter = vals.find(
"type");
1059 if ( iter != vals.end() ) {
1060 std::string type = iter->second;
1061 std::string tag = type +
"_Geant4_action";
1062 long result = PluginService::Create<long>(tag, &
m_detDesc, ptr, &vals);
1065 result = PluginService::Create<long>(tag, &
m_detDesc, ptr, &vals);
1067 except(
"LCDDConverter",
"Failed to locate plugin to interprete files of type"
1068 " \"" + tag +
"\" - no factory:" + type +
". " +
1072 result = *(
long*) result;
1074 except(
"LCDDConverter",
"Failed to invoke the plugin " + tag +
" of type " + type);
1076 printout(INFO,
"",
"+++ Executed Successfully Detector setup module %s.", type.c_str());
1079 printout(INFO,
"",
"+++ FAILED to execute Detector setup module %s.", nam.c_str());
1102 printout(WARNING,
"LCDDConverter",
"+++ No Detector header information availible from the geometry description.");
1105 template <
typename O,
typename C,
typename F>
void handle(
const O* o,
const C& c, F pmf) {
1106 for (
typename C::const_iterator i = c.begin(); i != c.end(); ++i) {
1107 std::string n = (*i)->GetName();
1112 template <
typename O,
typename C,
typename F>
void handleMap(
const O* o,
const C& c, F pmf) {
1113 for (
typename C::const_iterator i = c.begin(); i != c.end(); ++i)
1114 (o->*pmf)((*i).first, (*i).second);
1117 template <
typename O,
typename C,
typename F>
void handleRMap(
const O* o,
const C& c, F pmf) {
1118 for (
typename C::const_reverse_iterator i = c.rbegin(); i != c.rend(); ++i)
1119 handle(o, (*i).second, pmf);
1126 throw std::runtime_error(
"Attempt to call createGDML with an invalid geometry!");
1132 printout(ALWAYS,
"LCDDConverter",
"++ ==> Converting in memory detector description to GDML format...");
1138 "http://service-spi.web.cern.ch/service-spi/app/releases/GDML/schema/gdml.xsd");
1156 for(Detector::HandleMap::const_iterator i=mat.begin(); i!=mat.end(); ++i)
1157 geo.
materials.insert(
dynamic_cast<TGeoMedium*
>((*i).second.ptr()));
1162 printout(ALWAYS,
"LCDDConverter",
"++ Handled %ld materials.",geo.
materials.size());
1165 printout(ALWAYS,
"LCDDConverter",
"++ Collected %ld volumes.",geo.
volumes.size());
1168 printout(ALWAYS,
"LCDDConverter",
"++ Handled %ld solids.",geo.
solids.size());
1171 printout(ALWAYS,
"LCDDConverter",
"++ Handled %ld volumes.",geo.
volumes.size());
1181 throw std::runtime_error(
"Attempt to call createDetector with an invalid geometry!");
1187 printout(ALWAYS,
"LCDDConverter",
"++ ==> Dump visualisation attributes "
1188 "from in memory detector description...");
1198 printout(ALWAYS,
"LCDDConverter",
"++ Handled %ld volumes.",geo.
volumes.size());
1206 throw std::runtime_error(
"Attempt to call createDetector with an invalid geometry!");
1220 "http://www.lcsim.org/schemas/description/1.0/description.xsd");
1242 for (Detector::HandleMap::const_iterator i = fld.begin(); i != fld.end(); ++i)
1243 geo.
fields.insert((*i).second);
1245 printout(ALWAYS,
"LCDDConverter",
"++ ==> Converting in memory detector description to Detector format...");
1249 printout(ALWAYS,
"LCDDConverter",
"++ Handled %ld materials.",geo.
materials.size());
1252 printout(ALWAYS,
"LCDDConverter",
"++ Collected %ld volumes.",geo.
volumes.size());
1255 printout(ALWAYS,
"LCDDConverter",
"++ Handled %ld solids.",geo.
solids.size());
1258 printout(ALWAYS,
"LCDDConverter",
"++ Handled %ld visualization attributes.",geo.
vis.size());
1261 printout(ALWAYS,
"LCDDConverter",
"++ Handled %ld sensitive detectors.",geo.
sensitives.size());
1264 printout(ALWAYS,
"LCDDConverter",
"++ Handled %ld limit sets.",geo.
limits.size());
1267 printout(ALWAYS,
"LCDDConverter",
"++ Handled %ld regions.",geo.
regions.size());
1270 printout(ALWAYS,
"LCDDConverter",
"++ Handled %ld volumes.",geo.
volumes.size());
1273 printout(ALWAYS,
"LCDDConverter",
"++ Handled %ld fields.",geo.
fields.size());
1286 : doc(0), doc_root(0), doc_header(0), doc_idDict(0), doc_detectors(0), doc_limits(0),
1287 doc_regions(0), doc_display(0), doc_gdml(0), doc_fields(0), doc_define(0),
1288 doc_materials(0), doc_solids(0), doc_structure(0), doc_setup(0)
1292 static long dump_output(
xml_doc_t doc,
int argc,
char** argv) {
1294 return docH.
output(doc, argc > 0 ? argv[0] :
"");
1302 static long create_description(
Detector& description,
int argc,
char** argv) {
1304 return dump_output(wr.createDetector(description.
world()), argc, argv);
1307 static long create_vis(
Detector& description,
int argc,
char** argv) {
1309 return dump_output(wr.createVis(description.
world()), argc, argv);
1312 static long create_visASCII(
Detector& description,
int ,
char** argv) {
1314 wr.createVis(description.
world());
1316 std::map<std::string, xml_comp_t> vis_map;
1320 const char* sep =
";";
1321 std::ofstream os(argv[0]);
1325 auto iter = vis_map.find(ref.refStr());
1326 if ( iter != vis_map.end() ) {
1329 os <<
"vol:" << vol.nameStr() << sep <<
"vis:" << vis.nameStr() << sep
1330 <<
"visible:" << vis.visible() << sep <<
"r:"
1331 << col.R() << sep <<
"g:" << col.G() << sep <<
"b:" << col.B() << sep
1332 <<
"alpha:" << col.alpha() << sep <<
"line_style:"
1333 << vis.attr < std::string > (
_U(line_style)) << sep
1334 <<
"drawing_style:" << vis.attr < std::string> (
_U(drawing_style)) << sep
1335 <<
"show_daughters:" << vis.show_daughters() << sep << std::endl;