Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
TGDMLWrite.cxx
Go to the documentation of this file.
1 // @(#)root/gdml:$Id$
2 // Author: Anton Pytel 15/9/2011
3 
4 /*************************************************************************
5  * Copyright (C) 1995-2011, Rene Brun and Fons Rademakers. *
6  * All rights reserved. *
7  * *
8  * For the licensing terms see $ROOTSYS/LICENSE. *
9  * For the list of contributors see $ROOTSYS/README/CREDITS. *
10  *************************************************************************/
11 
12 /** \class TGDMLWrite
13 \ingroup Geometry_gdml
14 
15  This class contains implementation of converting ROOT's gGeoManager
16 geometry to GDML file. gGeoManager is the instance of TGeoManager class
17 containing tree of geometries creating resulting geometry. GDML is xml
18 based format of file mirroring the tree of geometries according to GDML
19 schema rules. For more information about GDML see http://gdml.web.cern.ch.
20 Each object in ROOT is represented by xml tag (=xml node/element) in GDML.
21 
22  This class is not needed to be instanciated. It should always be called
23 by gGeoManager->Export("xyz.gdml") method. Export is driven by extenstion
24 that is why ".gdml" is important in resulting name.
25 
26  Whenever a new ROOT geometry object is implemented or there is a change
27 in GDML schema this class is needed to be updated to ensure proper mapping
28 between ROOT objects and GDML elements.
29 
30  Current status of mapping ROOT -> GDML is implemented in method called
31 TGDMLWrite::ChooseObject and it contains following "map":
32 
33 #### Solids:
34 
35 ~~~
36 TGeoBBox -> <box ... >
37 TGeoParaboloid -> <paraboloid ...>
38 TGeoSphere -> <sphere ...>
39 TGeoArb8 -> <arb8 ...>
40 TGeoConeSeg -> <cone ...>
41 TGeoCone -> <cone ...>
42 TGeoPara -> <para ...>
43 TGeoTrap -> <trap ...> or
44 - -> <arb8 ...>
45 TGeoGtra -> <twistedtrap ...> or
46 - -> <trap ...> or
47 - -> <arb8 ...>
48 TGeoTrd1 -> <trd ...>
49 TGeoTrd2 -> <trd ...>
50 TGeoTubeSeg -> <tube ...>
51 TGeoCtub -> <cutTube ...>
52 TGeoTube -> <tube ...>
53 TGeoPcon -> <polycone ...>
54 TGeoTorus -> <torus ...>
55 TGeoPgon -> <polyhedra ...>
56 TGeoEltu -> <eltube ...>
57 TGeoHype -> <hype ...>
58 TGeoXtru -> <xtru ...>
59 TGeoCompositeShape -> <union ...> or
60 - -> <subtraction ...> or
61 - -> <intersection ...>
62 
63 Special cases of solids:
64 TGeoScaledShape -> <elcone ...> if scaled TGeoCone or
65 - -> element without scale
66 TGeoCompositeShape -> <ellipsoid ...>
67 - intersection of:
68 - scaled TGeoSphere and TGeoBBox
69 ~~~
70 
71 #### Materials:
72 
73 ~~~
74 TGeoIsotope -> <isotope ...>
75 TGeoElement -> <element ...>
76 TGeoMaterial -> <material ...>
77 TGeoMixture -> <material ...>
78 ~~~
79 
80 #### Structure
81 
82 ~~~
83 TGeoVolume -> <volume ...> or
84 - -> <assembly ...>
85 TGeoNode -> <physvol ...>
86 TGeoPatternFinder -> <divisionvol ...>
87 ~~~
88 
89 There are options that can be set to change resulting document
90 
91 ##### Options:
92 
93 ~~~
94 g - is set by default in gGeoManager, this option ensures compatibility
95 - with Geant4. It means:
96 - -> atomic number of material will be changed if <1 to 1
97 - -> if polycone is set badly it will try to export it correctly
98 - -> if widht * ndiv + offset is more then width of object being divided
99 - (in divisions) then it will be rounded so it will not exceed or
100 - if kPhi divsion then it will keep range of offset in -360 -> 0
101 f - if this option is set then names of volumes and solids will have
102 - pointer as a suffix to ensure uniqness of names
103 n - if this option is set then names will not have suffix, but uniqness is
104 - of names is not secured
105 - - if none of this two options (f,n) is set then default behaviour is so
106 - that incremental suffix is added to the names.
107 - (eg. TGeoBBox_0x1, TGeoBBox_0x2 ...)
108 ~~~
109 
110 #### USAGE:
111 
112 ~~~
113 gGeoManager->Export("output.gdml");
114 gGeoManager->Export("output.gdml","","vg"); //the same as previous just
115  //options are set explicitly
116 gGeoManager->Export("output.gdml","","vgf");
117 gGeoManager->Export("output.gdml","","gn");
118 gGeoManager->Export("output.gdml","","f");
119 ...
120 ~~~
121 
122 #### Note:
123  Options discussed above are used only for TGDMLWrite class. There are
124 other options in the TGeoManager::Export(...) method that can be used.
125 See that function for details.
126 
127 */
128 
129 #include "TGDMLWrite.h"
130 
131 #include "TGeoManager.h"
132 #include "TGeoMaterial.h"
133 #include "TGeoMatrix.h"
134 #include "TXMLEngine.h"
135 #include "TGeoVolume.h"
136 #include "TGeoBBox.h"
137 #include "TGeoParaboloid.h"
138 #include "TGeoArb8.h"
139 #include "TGeoTube.h"
140 #include "TGeoCone.h"
141 #include "TGeoTrd1.h"
142 #include "TGeoTrd2.h"
143 #include "TGeoPcon.h"
144 #include "TGeoPgon.h"
145 #include "TGeoSphere.h"
146 #include "TGeoTorus.h"
147 #include "TGeoPara.h"
148 #include "TGeoHype.h"
149 #include "TGeoEltu.h"
150 #include "TGeoXtru.h"
151 #include "TGeoScaledShape.h"
152 #include "TGeoVolume.h"
153 #include "TROOT.h"
154 #include "TMath.h"
155 #include "TGeoBoolNode.h"
156 #include "TGeoMedium.h"
157 #include "TGeoElement.h"
158 #include "TGeoShape.h"
159 #include "TGeoCompositeShape.h"
160 #include "TGeoOpticalSurface.h"
161 #include <stdlib.h>
162 #include <string>
163 #include <map>
164 #include <set>
165 #include <ctime>
166 #include <sstream>
167 
168 ClassImp(TGDMLWrite);
169 
170 TGDMLWrite *TGDMLWrite::fgGDMLWrite = 0;
171 
172 namespace {
173  struct MaterialExtractor {
174  std::set<TGeoMaterial*> materials;
175  void operator() (const TGeoVolume* v) {
176  materials.insert(v->GetMaterial());
177  for(Int_t i=0; i<v->GetNdaughters(); ++i)
178  (*this)(v->GetNode(i)->GetVolume());
179  }
180  };
181 }
182 
183 ////////////////////////////////////////////////////////////////////////////////
184 /// Default constructor.
185 
186 TGDMLWrite::TGDMLWrite()
187  : TObject(),
188  fIsotopeList(0),
189  fElementList(0),
190  fAccPatt(0),
191  fRejShape(0),
192  fNameList(0),
193  fgNamingSpeed(0),
194  fgG4Compatibility(0),
195  fGdmlFile(0),
196  fTopVolumeName(0),
197  fGdmlE(0),
198  fDefineNode(0),
199  fMaterialsNode(0),
200  fSolidsNode(0),
201  fStructureNode(0),
202  fVolCnt(0),
203  fPhysVolCnt(0),
204  fActNameErr(0),
205  fSolCnt(0),
206  fFltPrecision(17) // %.17g
207 {
208  if (fgGDMLWrite) delete fgGDMLWrite;
209  fgGDMLWrite = this;
210 }
211 
212 ////////////////////////////////////////////////////////////////////////////////
213 /// Destructor.
214 
215 TGDMLWrite::~TGDMLWrite()
216 {
217  delete fIsotopeList;
218  delete fElementList;
219  delete fAccPatt;
220  delete fRejShape;
221  delete fNameList;
222 
223  fgGDMLWrite = 0;
224 }
225 
226 ////////////////////////////////////////////////////////////////////////////////
227 /// Set convention of naming solids and volumes
228 
229 void TGDMLWrite::SetNamingSpeed(ENamingType naming)
230 {
231  fgNamingSpeed = naming;
232 }
233 
234 ////////////////////////////////////////////////////////////////////////////////
235 //wrapper of all main methods for extraction
236 void TGDMLWrite::WriteGDMLfile(TGeoManager * geomanager, const char* filename, TString option)
237 {
238  TList* materials = geomanager->GetListOfMaterials();
239  TGeoNode* node = geomanager->GetTopNode();
240  if ( !node ) {
241  Info("WriteGDMLfile", "Top volume does not exist!");
242  return;
243  }
244  fTopVolumeName = "";
245  WriteGDMLfile(geomanager, node, materials, filename, option);
246 }
247 
248 ////////////////////////////////////////////////////////////////////////////////
249 // Wrapper to only selectively write one branch of the volume hierarchy to file
250 void TGDMLWrite::WriteGDMLfile(TGeoManager * geomanager, TGeoNode* node, const char* filename, TString option)
251 {
252  TGeoVolume* volume = node->GetVolume();
253  TList materials, volumes, nodes;
254  MaterialExtractor extract;
255  if ( !volume ) {
256  Info("WriteGDMLfile", "Invalid Volume reference to extract GDML information!");
257  return;
258  }
259  extract(volume);
260  for(TGeoMaterial* m : extract.materials)
261  materials.Add(m);
262  fTopVolumeName = volume->GetName();
263  fSurfaceList.clear();
264  fVolumeList.clear();
265  fNodeList.clear();
266  WriteGDMLfile(geomanager, node, &materials, filename, option);
267  materials.Clear("nodelete");
268  volumes.Clear("nodelete");
269  nodes.Clear("nodelete");
270 }
271 
272 ////////////////////////////////////////////////////////////////////////////////
273 /// Wrapper of all exporting methods
274 /// Creates blank GDML file and fills it with gGeoManager structure converted
275 /// to GDML structure of xml nodes
276 
277 void TGDMLWrite::WriteGDMLfile(TGeoManager * geomanager,
278  TGeoNode* node,
279  TList* materialsLst,
280  const char* filename,
281  TString option)
282 {
283  //option processing
284  option.ToLower();
285  if (option.Contains("g")) {
286  SetG4Compatibility(kTRUE);
287  Info("WriteGDMLfile", "Geant4 compatibility mode set");
288  } else {
289  SetG4Compatibility(kFALSE);
290  }
291  if (option.Contains("f")) {
292  SetNamingSpeed(kfastButUglySufix);
293  Info("WriteGDMLfile", "Fast naming convention with pointer suffix set");
294  } else if (option.Contains("n")) {
295  SetNamingSpeed(kwithoutSufixNotUniq);
296  Info("WriteGDMLfile", "Naming without prefix set - be careful uniqness of name is not ensured");
297  } else {
298  SetNamingSpeed(kelegantButSlow);
299  Info("WriteGDMLfile", "Potentially slow with incremental suffix naming convention set");
300  }
301 
302  //local variables
303  Int_t outputLayout = 1;
304  const char * krootNodeName = "gdml";
305  const char * knsRefGeneral = "http://www.w3.org/2001/XMLSchema-instance";
306  const char * knsNameGeneral = "xsi";
307  const char * knsRefGdml = "http://service-spi.web.cern.ch/service-spi/app/releases/GDML/schema/gdml.xsd";
308  const char * knsNameGdml = "xsi:noNamespaceSchemaLocation";
309 
310  // First create engine
311  fGdmlE = new TXMLEngine;
312  fGdmlE->SetSkipComments(kTRUE);
313 
314  //create blank GDML file
315  fGdmlFile = fGdmlE->NewDoc();
316 
317  //create root node and add it to blank GDML file
318  XMLNodePointer_t rootNode = fGdmlE->NewChild(nullptr, nullptr, krootNodeName, nullptr);
319  fGdmlE->DocSetRootElement(fGdmlFile, rootNode);
320 
321  //add namespaces to root node
322  fGdmlE->NewNS(rootNode, knsRefGeneral, knsNameGeneral);
323  fGdmlE->NewAttr(rootNode, nullptr, knsNameGdml, knsRefGdml);
324 
325  //initialize general lists and <define>, <solids>, <structure> nodes
326  fIsotopeList = new StructLst;
327  fElementList = new StructLst;
328 
329  fNameList = new NameLst;
330 
331  fDefineNode = fGdmlE->NewChild(nullptr, nullptr, "define", nullptr);
332  fSolidsNode = fGdmlE->NewChild(nullptr, nullptr, "solids", nullptr);
333  fStructureNode = fGdmlE->NewChild(nullptr, nullptr, "structure", nullptr);
334  //========================
335 
336  //initialize list of accepted patterns for divisions (in ExtractVolumes)
337  fAccPatt = new StructLst;
338  fAccPatt->fLst["TGeoPatternX"] = kTRUE;
339  fAccPatt->fLst["TGeoPatternY"] = kTRUE;
340  fAccPatt->fLst["TGeoPatternZ"] = kTRUE;
341  fAccPatt->fLst["TGeoPatternCylR"] = kTRUE;
342  fAccPatt->fLst["TGeoPatternCylPhi"] = kTRUE;
343  //========================
344 
345  //initialize list of rejected shapes for divisions (in ExtractVolumes)
346  fRejShape = new StructLst;
347  //this shapes are rejected because, it is not possible to divide trd2
348  //in Y axis and while only trd2 object is imported from GDML
349  //it causes a problem when TGeoTrd1 is divided in Y axis
350  fRejShape->fLst["TGeoTrd1"] = kTRUE;
351  fRejShape->fLst["TGeoTrd2"] = kTRUE;
352  //=========================
353 
354  //Initialize global counters
355  fActNameErr = 0;
356  fVolCnt = 0;
357  fPhysVolCnt = 0;
358  fSolCnt = 0;
359 
360  //calling main extraction functions (with measuring time)
361  time_t startT, endT;
362  startT = time(nullptr);
363  ExtractMatrices(geomanager->GetListOfGDMLMatrices());
364  ExtractConstants(geomanager);
365  fMaterialsNode = ExtractMaterials(materialsLst);
366 
367  Info("WriteGDMLfile", "Extracting volumes");
368  ExtractVolumes(node);
369  Info("WriteGDMLfile", "%i solids added", fSolCnt);
370  Info("WriteGDMLfile", "%i volumes added", fVolCnt);
371  Info("WriteGDMLfile", "%i physvolumes added", fPhysVolCnt);
372  ExtractSkinSurfaces(geomanager->GetListOfSkinSurfaces());
373  ExtractBorderSurfaces(geomanager->GetListOfBorderSurfaces());
374  ExtractOpticalSurfaces(geomanager->GetListOfOpticalSurfaces());
375  endT = time(nullptr);
376  //<gdml>
377  fGdmlE->AddChild(rootNode, fDefineNode); // <define>...</define>
378  fGdmlE->AddChild(rootNode, fMaterialsNode); // <materials>...</materials>
379  fGdmlE->AddChild(rootNode, fSolidsNode); // <solids>...</solids>
380  fGdmlE->AddChild(rootNode, fStructureNode); // <structure>...</structure>
381  fGdmlE->AddChild(rootNode, CreateSetupN(fTopVolumeName.Data())); // <setup>...</setup>
382  //</gdml>
383  Double_t tdiffI = difftime(endT, startT);
384  TString tdiffS = (tdiffI == 0 ? TString("< 1 s") : TString::Format("%.0lf s", tdiffI));
385  Info("WriteGDMLfile", "Exporting time: %s", tdiffS.Data());
386  //=========================
387 
388  //Saving document
389  fGdmlE->SaveDoc(fGdmlFile, filename, outputLayout);
390  Info("WriteGDMLfile", "File %s saved", filename);
391  //cleaning
392  fGdmlE->FreeDoc(fGdmlFile);
393  //unset processing bits:
394  UnsetTemporaryBits(geomanager);
395  delete fGdmlE;
396 }
397 
398 ////////////////////////////////////////////////////////////////////////////////
399 /// Method exporting GDML matrices
400 
401 void TGDMLWrite::ExtractMatrices(TObjArray* matrixList)
402 {
403  if (!matrixList->GetEntriesFast()) return;
404  XMLNodePointer_t matrixN;
405  TIter next(matrixList);
406  TGDMLMatrix *matrix;
407  while ((matrix = (TGDMLMatrix*)next())) {
408  matrixN = CreateMatrixN(matrix);
409  fGdmlE->AddChild(fDefineNode, matrixN);
410  }
411 }
412 
413 ////////////////////////////////////////////////////////////////////////////////
414 /// Method exporting GDML matrices
415 
416 void TGDMLWrite::ExtractConstants(TGeoManager *geom)
417 {
418  if (!geom->GetNproperties()) return;
419  XMLNodePointer_t constantN;
420  TString property;
421  Double_t value;
422  for (Int_t i = 0; i < geom->GetNproperties(); ++i) {
423  value = geom->GetProperty(i, property);
424  constantN = CreateConstantN(property.Data(), value);
425  fGdmlE->AddChild(fDefineNode, constantN);
426  }
427 }
428 
429 ////////////////////////////////////////////////////////////////////////////////
430 /// Method exporting optical surfaces
431 
432 void TGDMLWrite::ExtractOpticalSurfaces(TObjArray *surfaces)
433 {
434  if (!surfaces->GetEntriesFast()) return;
435  XMLNodePointer_t surfaceN;
436  TIter next(surfaces);
437  TGeoOpticalSurface *surf;
438  while ((surf = (TGeoOpticalSurface*)next())) {
439  if ( fSurfaceList.find(surf) == fSurfaceList.end() ) continue;
440  surfaceN = CreateOpticalSurfaceN(surf);
441  fGdmlE->AddChild(fSolidsNode, surfaceN);
442  // Info("ExtractSkinSurfaces", "Extracted optical surface: %s",surf->GetName());
443  }
444 }
445 
446 ////////////////////////////////////////////////////////////////////////////////
447 /// Method exporting skin surfaces
448 
449 void TGDMLWrite::ExtractSkinSurfaces(TObjArray *surfaces)
450 {
451  if (!surfaces->GetEntriesFast()) return;
452  XMLNodePointer_t surfaceN;
453  TIter next(surfaces);
454  TGeoSkinSurface *surf;
455  while ((surf = (TGeoSkinSurface*)next())) {
456  if ( fVolumeList.find(surf->GetVolume()) == fVolumeList.end() ) continue;
457  surfaceN = CreateSkinSurfaceN(surf);
458  fGdmlE->AddChild(fStructureNode, surfaceN);
459  fSurfaceList.insert(surf->GetSurface());
460  // Info("ExtractSkinSurfaces", "Extracted skin surface: %s",surf->GetName());
461  }
462 }
463 
464 ////////////////////////////////////////////////////////////////////////////////
465 /// Method exporting border surfaces
466 
467 void TGDMLWrite::ExtractBorderSurfaces(TObjArray *surfaces)
468 {
469  if (!surfaces->GetEntriesFast()) return;
470  XMLNodePointer_t surfaceN;
471  TIter next(surfaces);
472  TGeoBorderSurface *surf;
473  while ((surf = (TGeoBorderSurface*)next())) {
474  auto ia = fNodeList.find(surf->GetNode1());
475  auto ib = fNodeList.find(surf->GetNode2());
476  if ( ia == fNodeList.end() && ib == fNodeList.end() ) {
477  continue;
478  }
479  else if ( ia == fNodeList.end() && ib != fNodeList.end() ) {
480  Warning("ExtractBorderSurfaces", "Inconsistent border surface extraction %s: Node %s"
481  " is not part of GDML!",surf->GetName(), surf->GetNode1()->GetName());
482  continue;
483  }
484  else if ( ia != fNodeList.end() && ib == fNodeList.end() ) {
485  Warning("ExtractBorderSurfaces", "Inconsistent border surface extraction %s: Node %s"
486  " is not part of GDML!",surf->GetName(), surf->GetNode2()->GetName());
487  continue;
488  }
489  surfaceN = CreateBorderSurfaceN(surf);
490  fGdmlE->AddChild(fStructureNode, surfaceN);
491  fSurfaceList.insert(surf->GetSurface());
492  // Info("ExtractBorderSurfaces", "Extracted border surface: %s",surf->GetName());
493  }
494 }
495 
496 ////////////////////////////////////////////////////////////////////////////////
497 /// Method exporting materials
498 
499 XMLNodePointer_t TGDMLWrite::ExtractMaterials(TList* materialsLst)
500 {
501  Info("ExtractMaterials", "Extracting materials");
502  //crate main <materials> node
503  XMLNodePointer_t materialsN = fGdmlE->NewChild(nullptr, nullptr, "materials", nullptr);
504  Int_t matcnt = 0;
505 
506  //go through materials - iterator and object declaration
507  TIter next(materialsLst);
508  TGeoMaterial *lmaterial;
509 
510  while ((lmaterial = (TGeoMaterial *)next())) {
511  //generate uniq name
512  TString lname = GenName(lmaterial->GetName(), TString::Format("%p", lmaterial));
513 
514  if (lmaterial->IsMixture()) {
515  TGeoMixture *lmixture = (TGeoMixture *)lmaterial;
516  XMLNodePointer_t mixtureN = CreateMixtureN(lmixture, materialsN, lname);
517  fGdmlE->AddChild(materialsN, mixtureN);
518  } else {
519  XMLNodePointer_t materialN = CreateMaterialN(lmaterial, lname);
520  fGdmlE->AddChild(materialsN, materialN);
521  }
522  matcnt++;
523  }
524  Info("ExtractMaterials", "%i materials added", matcnt);
525  return materialsN;
526 }
527 
528 ////////////////////////////////////////////////////////////////////////////////
529 /// Method creating solid to xml file and returning its name
530 
531 TString TGDMLWrite::ExtractSolid(TGeoShape* volShape)
532 {
533  XMLNodePointer_t solidN;
534  TString solname = "";
535  solidN = ChooseObject(volShape); //volume->GetShape()
536  fGdmlE->AddChild(fSolidsNode, solidN);
537  if (solidN != nullptr) fSolCnt++;
538  solname = fNameList->fLst[TString::Format("%p", volShape)];
539  if (solname.Contains("missing_")) {
540  solname = "-1";
541  }
542  return solname;
543 }
544 
545 
546 ////////////////////////////////////////////////////////////////////////////////
547 /// Method extracting geometry structure recursively
548 
549 void TGDMLWrite::ExtractVolumes(TGeoNode* node)
550 {
551  XMLNodePointer_t volumeN, childN;
552  TGeoVolume * volume = node->GetVolume();
553  TString volname, matname, solname, pattClsName, nodeVolNameBak;
554  TGeoPatternFinder *pattFinder = 0;
555  Bool_t isPattern = kFALSE;
556  const TString fltPrecision = TString::Format("%%.%dg", fFltPrecision);
557 
558  fNodeList.insert(node);
559  fVolumeList.insert(volume);
560  //create the name for volume/assembly
561  if (volume->IsTopVolume()) {
562  //not needed a special function for generating name
563  volname = volume->GetName();
564  fTopVolumeName = volname;
565  //register name to the pointer
566  fNameList->fLst[TString::Format("%p", volume)] = volname;
567  } else {
568  volname = GenName(volume->GetName(), TString::Format("%p", volume));
569  }
570 
571  //start to create main volume/assembly node
572  if (volume->IsAssembly()) {
573  volumeN = StartAssemblyN(volname);
574  } else {
575  //get reference material and add solid to <solids> + get name
576  matname = fNameList->fLst[TString::Format("%p", volume->GetMaterial())];
577  solname = ExtractSolid(volume->GetShape());
578  //If solid is not supported or corrupted
579  if (solname == "-1") {
580  Info("ExtractVolumes", "ERROR! %s volume was not added, because solid is either not supported or corrupted",
581  volname.Data());
582  //set volume as missing volume
583  fNameList->fLst[TString::Format("%p", volume)] = "missing_" + volname;
584  return;
585  }
586  volumeN = StartVolumeN(volname, solname, matname);
587 
588  //divisionvol can't be in assembly
589  pattFinder = volume->GetFinder();
590  //if found pattern
591  if (pattFinder) {
592  pattClsName = TString::Format("%s", pattFinder->ClassName());
593  TString shapeCls = TString::Format("%s", volume->GetShape()->ClassName());
594  //if pattern in accepted pattern list and not in shape rejected list
595  if ((fAccPatt->fLst[pattClsName] == kTRUE) &&
596  (fRejShape->fLst[shapeCls] != kTRUE)) {
597  isPattern = kTRUE;
598  }
599  }
600  }
601  //get all nodes in volume
602  TObjArray *nodeLst = volume->GetNodes();
603  TIter next(nodeLst);
604  TGeoNode *geoNode;
605  Int_t nCnt = 0;
606  //loop through all nodes
607  while ((geoNode = (TGeoNode *) next())) {
608  //get volume of current node and if not processed then process it
609  TGeoVolume * subvol = geoNode->GetVolume();
610  fNodeList.insert(geoNode);
611  if (subvol->TestAttBit(fgkProcBitVol) == kFALSE) {
612  subvol->SetAttBit(fgkProcBitVol);
613  ExtractVolumes(geoNode);
614  }
615 
616  //volume of this node has to exist because it was processed recursively
617  TString nodevolname = fNameList->fLst[TString::Format("%p", geoNode->GetVolume())];
618  if (nodevolname.Contains("missing_")) {
619  continue;
620  }
621  if (nCnt == 0) { //save name of the first node for divisionvol
622  nodeVolNameBak = nodevolname;
623  }
624 
625  if (isPattern == kFALSE) {
626  //create name for node
627  TString nodename, posname, rotname;
628  nodename = GenName(geoNode->GetName(), TString::Format("%p", geoNode));
629  nodename = nodename + "in" + volname;
630 
631  //create name for position and clear rotation
632  posname = nodename + "pos";
633  rotname = "";
634 
635  //position
636  const Double_t * pos = geoNode->GetMatrix()->GetTranslation();
637  Xyz nodPos;
638  nodPos.x = pos[0];
639  nodPos.y = pos[1];
640  nodPos.z = pos[2];
641  childN = CreatePositionN(posname.Data(), nodPos);
642  fGdmlE->AddChild(fDefineNode, childN); //adding node to <define> node
643  //Deal with reflection
644  XMLNodePointer_t scaleN = nullptr;
645  Double_t lx, ly, lz;
646  Double_t xangle = 0;
647  Double_t zangle = 0;
648  lx = geoNode->GetMatrix()->GetRotationMatrix()[0];
649  ly = geoNode->GetMatrix()->GetRotationMatrix()[4];
650  lz = geoNode->GetMatrix()->GetRotationMatrix()[8];
651  if (geoNode->GetMatrix()->IsReflection()
652  && TMath::Abs(lx) == 1 && TMath::Abs(ly) == 1 && TMath::Abs(lz) == 1) {
653  scaleN = fGdmlE->NewChild(nullptr, nullptr, "scale", nullptr);
654  fGdmlE->NewAttr(scaleN, nullptr, "name", (nodename + "scl").Data());
655  fGdmlE->NewAttr(scaleN, nullptr, "x", TString::Format(fltPrecision.Data(), lx));
656  fGdmlE->NewAttr(scaleN, nullptr, "y", TString::Format(fltPrecision.Data(), ly));
657  fGdmlE->NewAttr(scaleN, nullptr, "z", TString::Format(fltPrecision.Data(), lz));
658  //experimentally found out, that rotation should be updated like this
659  if (lx == -1) {
660  zangle = 180;
661  }
662  if (lz == -1) {
663  xangle = 180;
664  }
665  }
666 
667  //rotation
668  TGDMLWrite::Xyz lxyz = GetXYZangles(geoNode->GetMatrix()->GetRotationMatrix());
669  lxyz.x -= xangle;
670  lxyz.z -= zangle;
671  if ((lxyz.x != 0.0) || (lxyz.y != 0.0) || (lxyz.z != 0.0)) {
672  rotname = nodename + "rot";
673  childN = CreateRotationN(rotname.Data(), lxyz);
674  fGdmlE->AddChild(fDefineNode, childN); //adding node to <define> node
675  }
676 
677  //create physvol for main volume/assembly node
678  childN = CreatePhysVolN(geoNode->GetName(), geoNode->GetNumber(), nodevolname.Data(), posname.Data(), rotname.Data(), scaleN);
679  fGdmlE->AddChild(volumeN, childN);
680  }
681  nCnt++;
682  }
683  //create only one divisionvol node
684  if (isPattern && pattFinder) {
685  //retrieve attributes of division
686  Int_t ndiv, divaxis;
687  Double_t offset, width, xlo, xhi;
688  TString axis, unit;
689 
690  ndiv = pattFinder->GetNdiv();
691  width = pattFinder->GetStep();
692 
693  divaxis = pattFinder->GetDivAxis();
694  volume->GetShape()->GetAxisRange(divaxis, xlo, xhi);
695 
696  //compute relative start (not positional)
697  offset = pattFinder->GetStart() - xlo;
698  axis = GetPattAxis(divaxis, pattClsName, unit);
699 
700  //create division node
701  childN = CreateDivisionN(offset, width, ndiv, axis.Data(), unit.Data(), nodeVolNameBak.Data());
702  fGdmlE->AddChild(volumeN, childN);
703  }
704 
705  fVolCnt++;
706  //add volume/assembly node into the <structure> node
707  fGdmlE->AddChild(fStructureNode, volumeN);
708 
709 }
710 
711 ////////////////////////////////////////////////////////////////////////////////
712 /// Creates "atom" node for GDML
713 
714 XMLNodePointer_t TGDMLWrite::CreateAtomN(Double_t atom, const char * unit)
715 {
716  const TString fltPrecision = TString::Format("%%.%dg", fFltPrecision);
717  XMLNodePointer_t atomN = fGdmlE->NewChild(nullptr, nullptr, "atom", nullptr);
718  fGdmlE->NewAttr(atomN, nullptr, "unit", unit);
719  fGdmlE->NewAttr(atomN, nullptr, "value", TString::Format(fltPrecision.Data(), atom));
720  return atomN;
721 }
722 
723 ////////////////////////////////////////////////////////////////////////////////
724 /// Creates "D" density node for GDML
725 
726 XMLNodePointer_t TGDMLWrite::CreateDN(Double_t density, const char * unit)
727 {
728  const TString fltPrecision = TString::Format("%%.%dg", fFltPrecision);
729  XMLNodePointer_t densN = fGdmlE->NewChild(nullptr, nullptr, "D", nullptr);
730  fGdmlE->NewAttr(densN, nullptr, "unit", unit);
731  fGdmlE->NewAttr(densN, nullptr, "value", TString::Format(fltPrecision.Data(), density));
732  return densN;
733 }
734 
735 ////////////////////////////////////////////////////////////////////////////////
736 /// Creates "fraction" node for GDML
737 
738 XMLNodePointer_t TGDMLWrite::CreateFractionN(Double_t percentage, const char * refName)
739 {
740  const TString fltPrecision = TString::Format("%%.%dg", fFltPrecision);
741  XMLNodePointer_t fractN = fGdmlE->NewChild(nullptr, nullptr, "fraction", nullptr);
742  fGdmlE->NewAttr(fractN, nullptr, "n", TString::Format(fltPrecision.Data(), percentage));
743  fGdmlE->NewAttr(fractN, nullptr, "ref", refName);
744  return fractN;
745 }
746 
747 ////////////////////////////////////////////////////////////////////////////////
748 /// Creates "property" node for GDML
749 
750 XMLNodePointer_t TGDMLWrite::CreatePropertyN(TNamed const &property)
751 {
752  XMLNodePointer_t propertyN = fGdmlE->NewChild(nullptr, nullptr, "property", nullptr);
753  fGdmlE->NewAttr(propertyN, nullptr, "name", property.GetName());
754  fGdmlE->NewAttr(propertyN, nullptr, "ref", property.GetTitle());
755  return propertyN;
756 }
757 
758 ////////////////////////////////////////////////////////////////////////////////
759 /// Creates "isotope" node for GDML
760 
761 XMLNodePointer_t TGDMLWrite::CreateIsotopN(TGeoIsotope * isotope, const char * name)
762 {
763  XMLNodePointer_t mainN = fGdmlE->NewChild(nullptr, nullptr, "isotope", nullptr);
764  fGdmlE->NewAttr(mainN, nullptr, "name", name);
765  fGdmlE->NewAttr(mainN, nullptr, "N", TString::Format("%i", isotope->GetN()));
766  fGdmlE->NewAttr(mainN, nullptr, "Z", TString::Format("%i", isotope->GetZ()));
767  fGdmlE->AddChild(mainN, CreateAtomN(isotope->GetA()));
768  return mainN;
769 }
770 
771 ////////////////////////////////////////////////////////////////////////////////
772 /// Creates "element" node for GDML
773 ///element node and attribute
774 
775 XMLNodePointer_t TGDMLWrite::CreateElementN(TGeoElement * element, XMLNodePointer_t materials, const char * name)
776 {
777  XMLNodePointer_t mainN = fGdmlE->NewChild(nullptr, nullptr, "element", nullptr);
778  fGdmlE->NewAttr(mainN, nullptr, "name", name);
779  //local associative arrays for saving isotopes and their weight
780  //inside element
781  NameListF wPercentage;
782  NameListI wCounter;
783 
784  if (element->HasIsotopes()) {
785  Int_t nOfIso = element->GetNisotopes();
786  //go through isotopes
787  for (Int_t idx = 0; idx < nOfIso; idx++) {
788  TGeoIsotope *myIsotope = element->GetIsotope(idx);
789  if (!myIsotope) {
790  Fatal("CreateElementN", "Missing isotopes for element %s", element->GetName());
791  return mainN;
792  }
793 
794  //Get name of the Isotope (
795  TString lname = myIsotope->GetName();
796  //_iso suffix is added to avoid problems with same names
797  //for material, element and isotopes
798  lname = TString::Format("%s_iso", lname.Data());
799 
800  //cumulates abundance, in case 2 isotopes with same names
801  //within one element
802  wPercentage[lname] += element->GetRelativeAbundance(idx);
803  wCounter[lname]++;
804 
805  //check whether isotope name is not in list of isotopes
806  if (IsInList(fIsotopeList->fLst, lname)) {
807  continue;
808  }
809  //add isotope to list of isotopes and to main <materials> node
810  fIsotopeList->fLst[lname] = kTRUE;
811  XMLNodePointer_t isoNode = CreateIsotopN(myIsotope, lname);
812  fGdmlE->AddChild(materials, isoNode);
813  }
814  //loop through asoc array of isotopes
815  for (NameListI::iterator itr = wCounter.begin(); itr != wCounter.end(); ++itr) {
816  if (itr->second > 1) {
817  Info("CreateMixtureN", "WARNING! 2 equal isotopes in one element. Check: %s isotope of %s element",
818  itr->first.Data(), name);
819  }
820  //add fraction child to element with reference to isotope
821  fGdmlE->AddChild(mainN, CreateFractionN(wPercentage[itr->first], itr->first.Data()));
822  }
823  } else {
824  fGdmlE->NewAttr(mainN, nullptr, "formula", element->GetName());
825  Int_t valZ = element->Z();
826  // Z can't be <1 in Geant4 and Z is optional parameter
827  if (valZ >= 1) {
828  fGdmlE->NewAttr(mainN, nullptr, "Z", TString::Format("%i", valZ));
829  }
830  fGdmlE->AddChild(mainN, CreateAtomN(element->A()));
831  }
832  return mainN;
833 }
834 
835 ////////////////////////////////////////////////////////////////////////////////
836 /// Creates "material" node for GDML with references to other sub elements
837 
838 XMLNodePointer_t TGDMLWrite::CreateMixtureN(TGeoMixture * mixture, XMLNodePointer_t materials, TString mname)
839 {
840  XMLNodePointer_t mainN = fGdmlE->NewChild(nullptr, nullptr, "material", nullptr);
841  fGdmlE->NewAttr(mainN, nullptr, "name", mname);
842  fGdmlE->AddChild(mainN, CreateDN(mixture->GetDensity()));
843  //local associative arrays for saving elements and their weight
844  //inside mixture
845  NameListF wPercentage;
846  NameListI wCounter;
847 
848  Int_t nOfElm = mixture->GetNelements();
849  //go through elements
850  for (Int_t idx = 0; idx < nOfElm; idx++) {
851  TGeoElement *myElement = mixture->GetElement(idx);
852 
853  //Get name of the element
854  //NOTE: that for element - GetTitle() returns the "name" tag
855  //and GetName() returns "formula" tag (see createElementN)
856  TString lname = myElement->GetTitle();
857  //_elm suffix is added to avoid problems with same names
858  //for material and element
859  lname = TString::Format("%s_elm", lname.Data());
860 
861  //cumulates percentage, in case 2 elements with same names within one mixture
862  wPercentage[lname] += mixture->GetWmixt()[idx];
863  wCounter[lname]++;
864 
865  //check whether element name is not in list of elements already created
866  if (IsInList(fElementList->fLst, lname)) {
867  continue;
868  }
869 
870  //add element to list of elements and to main <materials> node
871  fElementList->fLst[lname] = kTRUE;
872  XMLNodePointer_t elmNode = CreateElementN(myElement, materials, lname);
873  fGdmlE->AddChild(materials, elmNode);
874  }
875  //loop through asoc array
876  for (NameListI::iterator itr = wCounter.begin(); itr != wCounter.end(); ++itr) {
877  if (itr->second > 1) {
878  Info("CreateMixtureN", "WARNING! 2 equal elements in one material. Check: %s element of %s material",
879  itr->first.Data(), mname.Data());
880  }
881  //add fraction child to material with reference to element
882  fGdmlE->AddChild(mainN, CreateFractionN(wPercentage[itr->first], itr->first.Data()));
883  }
884 
885  // Write properties
886  TList const &properties = mixture->GetProperties();
887  if (properties.GetSize()) {
888  TIter next(&properties);
889  TNamed *property;
890  while ((property = (TNamed*)next()))
891  fGdmlE->AddChild(mainN, CreatePropertyN(*property));
892  }
893  // Write CONST properties
894  TList const &const_properties = mixture->GetConstProperties();
895  if (const_properties.GetSize()) {
896  TIter next(&const_properties);
897  TNamed *property;
898  while ((property = (TNamed*)next()))
899  fGdmlE->AddChild(mainN, CreatePropertyN(*property));
900  }
901 
902  return mainN;
903 }
904 
905 ////////////////////////////////////////////////////////////////////////////////
906 /// Creates "material" node for GDML
907 
908 XMLNodePointer_t TGDMLWrite::CreateMaterialN(TGeoMaterial * material, TString mname)
909 {
910  XMLNodePointer_t mainN = fGdmlE->NewChild(nullptr, nullptr, "material", nullptr);
911  fGdmlE->NewAttr(mainN, nullptr, "name", mname);
912  Double_t valZ = material->GetZ();
913  //Z can't be zero in Geant4 so this is workaround for vacuum
914  TString tmpname = mname;
915  const TString fltPrecision = TString::Format("%%.%dg", fFltPrecision);
916  tmpname.ToLower();
917  if (valZ < 1) {
918  if (tmpname == "vacuum") {
919  valZ = 1;
920  } else {
921  if (fgG4Compatibility == kTRUE) {
922  Info("CreateMaterialN", "WARNING! value of Z in %s material can't be < 1 in Geant4, that is why it was changed to 1, please check it manually! ",
923  mname.Data());
924  valZ = 1;
925  } else {
926  Info("CreateMaterialN", "WARNING! value of Z in %s material can't be < 1 in Geant4", mname.Data());
927  }
928  }
929  }
930  fGdmlE->NewAttr(mainN, nullptr, "Z", TString::Format(fltPrecision.Data(), valZ)); //material->GetZ()));
931  fGdmlE->AddChild(mainN, CreateDN(material->GetDensity()));
932  fGdmlE->AddChild(mainN, CreateAtomN(material->GetA()));
933  // Create properties if any
934  TList const &properties = material->GetProperties();
935  if (properties.GetSize()) {
936  TIter next(&properties);
937  TNamed *property;
938  while ((property = (TNamed*)next()))
939  fGdmlE->AddChild(mainN, CreatePropertyN(*property));
940  }
941  // Write CONST properties
942  TList const &const_properties = material->GetConstProperties();
943  if (const_properties.GetSize()) {
944  TIter next(&const_properties);
945  TNamed *property;
946  while ((property = (TNamed*)next()))
947  fGdmlE->AddChild(mainN, CreatePropertyN(*property));
948  }
949  return mainN;
950 }
951 
952 ////////////////////////////////////////////////////////////////////////////////
953 /// Creates "box" node for GDML
954 
955 XMLNodePointer_t TGDMLWrite::CreateBoxN(TGeoBBox * geoShape)
956 {
957  XMLNodePointer_t mainN = fGdmlE->NewChild(nullptr, nullptr, "box", nullptr);
958  const TString fltPrecision = TString::Format("%%.%dg", fFltPrecision);
959  TString lname = GenName(geoShape->GetName(), TString::Format("%p", geoShape));
960  fGdmlE->NewAttr(mainN, nullptr, "name", lname);
961  if (IsNullParam(geoShape->GetDX(), "DX", lname) ||
962  IsNullParam(geoShape->GetDY(), "DY", lname) ||
963  IsNullParam(geoShape->GetDZ(), "DZ", lname)) {
964  return nullptr;
965  }
966  fGdmlE->NewAttr(mainN, nullptr, "x", TString::Format(fltPrecision.Data(), 2 * geoShape->GetDX()));
967  fGdmlE->NewAttr(mainN, nullptr, "y", TString::Format(fltPrecision.Data(), 2 * geoShape->GetDY()));
968  fGdmlE->NewAttr(mainN, nullptr, "z", TString::Format(fltPrecision.Data(), 2 * geoShape->GetDZ()));
969 
970  fGdmlE->NewAttr(mainN, nullptr, "lunit", "cm");
971  return mainN;
972 }
973 
974 ////////////////////////////////////////////////////////////////////////////////
975 /// Creates "paraboloid" node for GDML
976 
977 XMLNodePointer_t TGDMLWrite::CreateParaboloidN(TGeoParaboloid * geoShape)
978 {
979  XMLNodePointer_t mainN = fGdmlE->NewChild(nullptr, nullptr, "paraboloid", nullptr);
980  const TString fltPrecision = TString::Format("%%.%dg", fFltPrecision);
981  TString lname = GenName(geoShape->GetName(), TString::Format("%p", geoShape));
982  fGdmlE->NewAttr(mainN, nullptr, "name", lname);
983  if (IsNullParam(geoShape->GetRhi(), "Rhi", lname) ||
984  IsNullParam(geoShape->GetDz(), "Dz", lname)) {
985  return nullptr;
986  }
987  fGdmlE->NewAttr(mainN, nullptr, "rlo", TString::Format(fltPrecision.Data(), geoShape->GetRlo()));
988  fGdmlE->NewAttr(mainN, nullptr, "rhi", TString::Format(fltPrecision.Data(), geoShape->GetRhi()));
989  fGdmlE->NewAttr(mainN, nullptr, "dz", TString::Format(fltPrecision.Data(), geoShape->GetDz()));
990 
991  fGdmlE->NewAttr(mainN, nullptr, "lunit", "cm");
992  return mainN;
993 }
994 
995 ////////////////////////////////////////////////////////////////////////////////
996 /// Creates "sphere" node for GDML
997 
998 XMLNodePointer_t TGDMLWrite::CreateSphereN(TGeoSphere * geoShape)
999 {
1000  XMLNodePointer_t mainN = fGdmlE->NewChild(nullptr, nullptr, "sphere", nullptr);
1001  const TString fltPrecision = TString::Format("%%.%dg", fFltPrecision);
1002  TString lname = GenName(geoShape->GetName(), TString::Format("%p", geoShape));
1003  fGdmlE->NewAttr(mainN, nullptr, "name", lname);
1004  if (IsNullParam(geoShape->GetRmax(), "Rmax", lname)) {
1005  return nullptr;
1006  }
1007 
1008  fGdmlE->NewAttr(mainN, nullptr, "rmin", TString::Format(fltPrecision.Data(), geoShape->GetRmin()));
1009  fGdmlE->NewAttr(mainN, nullptr, "rmax", TString::Format(fltPrecision.Data(), geoShape->GetRmax()));
1010  fGdmlE->NewAttr(mainN, nullptr, "startphi", TString::Format(fltPrecision.Data(), geoShape->GetPhi1()));
1011  fGdmlE->NewAttr(mainN, nullptr, "deltaphi", TString::Format(fltPrecision.Data(), geoShape->GetPhi2() - geoShape->GetPhi1()));
1012  fGdmlE->NewAttr(mainN, nullptr, "starttheta", TString::Format(fltPrecision.Data(), geoShape->GetTheta1()));
1013  fGdmlE->NewAttr(mainN, nullptr, "deltatheta", TString::Format(fltPrecision.Data(), geoShape->GetTheta2() - geoShape->GetTheta1()));
1014 
1015  fGdmlE->NewAttr(mainN, nullptr, "aunit", "deg");
1016  fGdmlE->NewAttr(mainN, nullptr, "lunit", "cm");
1017  return mainN;
1018 }
1019 
1020 ////////////////////////////////////////////////////////////////////////////////
1021 /// Creates "arb8" node for GDML
1022 
1023 XMLNodePointer_t TGDMLWrite::CreateArb8N(TGeoArb8 * geoShape)
1024 {
1025  XMLNodePointer_t mainN = fGdmlE->NewChild(nullptr, nullptr, "arb8", nullptr);
1026  const TString fltPrecision = TString::Format("%%.%dg", fFltPrecision);
1027  TString lname = GenName(geoShape->GetName(), TString::Format("%p", geoShape));
1028  fGdmlE->NewAttr(mainN, nullptr, "name", lname);
1029  if (IsNullParam(geoShape->GetDz(), "Dz", lname)) {
1030  return nullptr;
1031  }
1032 
1033  fGdmlE->NewAttr(mainN, nullptr, "v1x", TString::Format(fltPrecision.Data(), geoShape->GetVertices()[0]));
1034  fGdmlE->NewAttr(mainN, nullptr, "v1y", TString::Format(fltPrecision.Data(), geoShape->GetVertices()[1]));
1035  fGdmlE->NewAttr(mainN, nullptr, "v2x", TString::Format(fltPrecision.Data(), geoShape->GetVertices()[2]));
1036  fGdmlE->NewAttr(mainN, nullptr, "v2y", TString::Format(fltPrecision.Data(), geoShape->GetVertices()[3]));
1037  fGdmlE->NewAttr(mainN, nullptr, "v3x", TString::Format(fltPrecision.Data(), geoShape->GetVertices()[4]));
1038  fGdmlE->NewAttr(mainN, nullptr, "v3y", TString::Format(fltPrecision.Data(), geoShape->GetVertices()[5]));
1039  fGdmlE->NewAttr(mainN, nullptr, "v4x", TString::Format(fltPrecision.Data(), geoShape->GetVertices()[6]));
1040  fGdmlE->NewAttr(mainN, nullptr, "v4y", TString::Format(fltPrecision.Data(), geoShape->GetVertices()[7]));
1041  fGdmlE->NewAttr(mainN, nullptr, "v5x", TString::Format(fltPrecision.Data(), geoShape->GetVertices()[8]));
1042  fGdmlE->NewAttr(mainN, nullptr, "v5y", TString::Format(fltPrecision.Data(), geoShape->GetVertices()[9]));
1043  fGdmlE->NewAttr(mainN, nullptr, "v6x", TString::Format(fltPrecision.Data(), geoShape->GetVertices()[10]));
1044  fGdmlE->NewAttr(mainN, nullptr, "v6y", TString::Format(fltPrecision.Data(), geoShape->GetVertices()[11]));
1045  fGdmlE->NewAttr(mainN, nullptr, "v7x", TString::Format(fltPrecision.Data(), geoShape->GetVertices()[12]));
1046  fGdmlE->NewAttr(mainN, nullptr, "v7y", TString::Format(fltPrecision.Data(), geoShape->GetVertices()[13]));
1047  fGdmlE->NewAttr(mainN, nullptr, "v8x", TString::Format(fltPrecision.Data(), geoShape->GetVertices()[14]));
1048  fGdmlE->NewAttr(mainN, nullptr, "v8y", TString::Format(fltPrecision.Data(), geoShape->GetVertices()[15]));
1049  fGdmlE->NewAttr(mainN, nullptr, "dz", TString::Format(fltPrecision.Data(), geoShape->GetDz()));
1050 
1051  fGdmlE->NewAttr(mainN, nullptr, "lunit", "cm");
1052  return mainN;
1053 }
1054 
1055 ////////////////////////////////////////////////////////////////////////////////
1056 /// Creates "cone" node for GDML from TGeoConeSeg object
1057 
1058 XMLNodePointer_t TGDMLWrite::CreateConeN(TGeoConeSeg * geoShape)
1059 {
1060  XMLNodePointer_t mainN = fGdmlE->NewChild(nullptr, nullptr, "cone", nullptr);
1061  const TString fltPrecision = TString::Format("%%.%dg", fFltPrecision);
1062  TString lname = GenName(geoShape->GetName(), TString::Format("%p", geoShape));
1063  fGdmlE->NewAttr(mainN, nullptr, "name", lname);
1064  if (IsNullParam(geoShape->GetDz(), "Dz", lname)) {
1065  return nullptr;
1066  }
1067 
1068  fGdmlE->NewAttr(mainN, nullptr, "z", TString::Format(fltPrecision.Data(), 2 * geoShape->GetDz()));
1069  fGdmlE->NewAttr(mainN, nullptr, "rmin1", TString::Format(fltPrecision.Data(), geoShape->GetRmin1()));
1070  fGdmlE->NewAttr(mainN, nullptr, "rmin2", TString::Format(fltPrecision.Data(), geoShape->GetRmin2()));
1071  fGdmlE->NewAttr(mainN, nullptr, "rmax1", TString::Format(fltPrecision.Data(), geoShape->GetRmax1()));
1072  fGdmlE->NewAttr(mainN, nullptr, "rmax2", TString::Format(fltPrecision.Data(), geoShape->GetRmax2()));
1073  fGdmlE->NewAttr(mainN, nullptr, "startphi", TString::Format(fltPrecision.Data(), geoShape->GetPhi1()));
1074  fGdmlE->NewAttr(mainN, nullptr, "deltaphi", TString::Format(fltPrecision.Data(), geoShape->GetPhi2() - geoShape->GetPhi1()));
1075 
1076  fGdmlE->NewAttr(mainN, nullptr, "aunit", "deg");
1077  fGdmlE->NewAttr(mainN, nullptr, "lunit", "cm");
1078  return mainN;
1079 }
1080 
1081 ////////////////////////////////////////////////////////////////////////////////
1082 /// Creates "cone" node for GDML from TGeoCone object
1083 
1084 XMLNodePointer_t TGDMLWrite::CreateConeN(TGeoCone * geoShape)
1085 {
1086  XMLNodePointer_t mainN = fGdmlE->NewChild(nullptr, nullptr, "cone", nullptr);
1087  const TString fltPrecision = TString::Format("%%.%dg", fFltPrecision);
1088  TString lname = GenName(geoShape->GetName(), TString::Format("%p", geoShape));
1089  fGdmlE->NewAttr(mainN, nullptr, "name", lname);
1090  if (IsNullParam(geoShape->GetDz(), "Dz", lname)) {
1091  return nullptr;
1092  }
1093 
1094  fGdmlE->NewAttr(mainN, nullptr, "z", TString::Format(fltPrecision.Data(), 2 * geoShape->GetDz()));
1095  fGdmlE->NewAttr(mainN, nullptr, "rmin1", TString::Format(fltPrecision.Data(), geoShape->GetRmin1()));
1096  fGdmlE->NewAttr(mainN, nullptr, "rmin2", TString::Format(fltPrecision.Data(), geoShape->GetRmin2()));
1097  fGdmlE->NewAttr(mainN, nullptr, "rmax1", TString::Format(fltPrecision.Data(), geoShape->GetRmax1()));
1098  fGdmlE->NewAttr(mainN, nullptr, "rmax2", TString::Format(fltPrecision.Data(), geoShape->GetRmax2()));
1099  fGdmlE->NewAttr(mainN, nullptr, "startphi", TString::Format("%i", 0));
1100  fGdmlE->NewAttr(mainN, nullptr, "deltaphi", TString::Format("%i", 360));
1101 
1102  fGdmlE->NewAttr(mainN, nullptr, "aunit", "deg");
1103  fGdmlE->NewAttr(mainN, nullptr, "lunit", "cm");
1104  return mainN;
1105 }
1106 
1107 ////////////////////////////////////////////////////////////////////////////////
1108 /// Creates "para" node for GDML
1109 
1110 XMLNodePointer_t TGDMLWrite::CreateParaN(TGeoPara * geoShape)
1111 {
1112  XMLNodePointer_t mainN = fGdmlE->NewChild(nullptr, nullptr, "para", nullptr);
1113  const TString fltPrecision = TString::Format("%%.%dg", fFltPrecision);
1114  fGdmlE->NewAttr(mainN, nullptr, "name", GenName(geoShape->GetName(), TString::Format("%p", geoShape)));
1115 
1116  fGdmlE->NewAttr(mainN, nullptr, "x", TString::Format(fltPrecision.Data(), 2 * geoShape->GetX()));
1117  fGdmlE->NewAttr(mainN, nullptr, "y", TString::Format(fltPrecision.Data(), 2 * geoShape->GetY()));
1118  fGdmlE->NewAttr(mainN, nullptr, "z", TString::Format(fltPrecision.Data(), 2 * geoShape->GetZ()));
1119  fGdmlE->NewAttr(mainN, nullptr, "alpha", TString::Format(fltPrecision.Data(), geoShape->GetAlpha()));
1120  fGdmlE->NewAttr(mainN, nullptr, "theta", TString::Format(fltPrecision.Data(), geoShape->GetTheta()));
1121  fGdmlE->NewAttr(mainN, nullptr, "phi", TString::Format(fltPrecision.Data(), geoShape->GetPhi()));
1122 
1123  fGdmlE->NewAttr(mainN, nullptr, "aunit", "deg");
1124  fGdmlE->NewAttr(mainN, nullptr, "lunit", "cm");
1125  return mainN;
1126 }
1127 
1128 ////////////////////////////////////////////////////////////////////////////////
1129 /// Creates "trap" node for GDML
1130 
1131 XMLNodePointer_t TGDMLWrite::CreateTrapN(TGeoTrap * geoShape)
1132 {
1133  XMLNodePointer_t mainN;
1134  const TString fltPrecision = TString::Format("%%.%dg", fFltPrecision);
1135 
1136  //if one base equals 0 create Arb8 instead of trap
1137  if ((geoShape->GetBl1() == 0 || geoShape->GetTl1() == 0 || geoShape->GetH1() == 0) ||
1138  (geoShape->GetBl2() == 0 || geoShape->GetTl2() == 0 || geoShape->GetH2() == 0)) {
1139  mainN = CreateArb8N(geoShape);
1140  return mainN;
1141  }
1142 
1143  //if is twisted then create arb8
1144  if (geoShape->IsTwisted()) {
1145  mainN = CreateArb8N((TGeoArb8 *) geoShape);
1146  return mainN;
1147  }
1148 
1149  mainN = fGdmlE->NewChild(nullptr, nullptr, "trap", nullptr);
1150  TString lname = GenName(geoShape->GetName(), TString::Format("%p", geoShape));
1151  fGdmlE->NewAttr(mainN, nullptr, "name", lname);
1152  if (IsNullParam(geoShape->GetDz(), "Dz", lname)) {
1153  return nullptr;
1154  }
1155 
1156  fGdmlE->NewAttr(mainN, nullptr, "z", TString::Format(fltPrecision.Data(), 2 * geoShape->GetDz()));
1157  fGdmlE->NewAttr(mainN, nullptr, "theta", TString::Format(fltPrecision.Data(), geoShape->GetTheta()));
1158  fGdmlE->NewAttr(mainN, nullptr, "phi", TString::Format(fltPrecision.Data(), geoShape->GetPhi()));
1159  fGdmlE->NewAttr(mainN, nullptr, "x1", TString::Format(fltPrecision.Data(), 2 * geoShape->GetBl1()));
1160  fGdmlE->NewAttr(mainN, nullptr, "x2", TString::Format(fltPrecision.Data(), 2 * geoShape->GetTl1()));
1161  fGdmlE->NewAttr(mainN, nullptr, "x3", TString::Format(fltPrecision.Data(), 2 * geoShape->GetBl2()));
1162  fGdmlE->NewAttr(mainN, nullptr, "x4", TString::Format(fltPrecision.Data(), 2 * geoShape->GetTl2()));
1163  fGdmlE->NewAttr(mainN, nullptr, "y1", TString::Format(fltPrecision.Data(), 2 * geoShape->GetH1()));
1164  fGdmlE->NewAttr(mainN, nullptr, "y2", TString::Format(fltPrecision.Data(), 2 * geoShape->GetH2()));
1165 
1166  fGdmlE->NewAttr(mainN, nullptr, "alpha1", TString::Format(fltPrecision.Data(), geoShape->GetAlpha1()));
1167  fGdmlE->NewAttr(mainN, nullptr, "alpha2", TString::Format(fltPrecision.Data(), geoShape->GetAlpha2()));
1168 
1169  fGdmlE->NewAttr(mainN, nullptr, "aunit", "deg");
1170  fGdmlE->NewAttr(mainN, nullptr, "lunit", "cm");
1171  return mainN;
1172 }
1173 
1174 ////////////////////////////////////////////////////////////////////////////////
1175 /// Creates "twistedtrap" node for GDML
1176 
1177 XMLNodePointer_t TGDMLWrite::CreateTwistedTrapN(TGeoGtra * geoShape)
1178 {
1179  XMLNodePointer_t mainN;
1180  const TString fltPrecision = TString::Format("%%.%dg", fFltPrecision);
1181 
1182  //if one base equals 0 create Arb8 instead of twisted trap
1183  if ((geoShape->GetBl1() == 0 && geoShape->GetTl1() == 0 && geoShape->GetH1() == 0) ||
1184  (geoShape->GetBl2() == 0 && geoShape->GetTl2() == 0 && geoShape->GetH2() == 0)) {
1185  mainN = CreateArb8N((TGeoArb8 *) geoShape);
1186  return mainN;
1187  }
1188 
1189  //if is twisted then create arb8
1190  if (geoShape->IsTwisted()) {
1191  mainN = CreateArb8N((TGeoArb8 *) geoShape);
1192  return mainN;
1193  }
1194 
1195  //if parameter twistAngle (PhiTwist) equals zero create trap node
1196  if (geoShape->GetTwistAngle() == 0) {
1197  mainN = CreateTrapN((TGeoTrap *) geoShape);
1198  return mainN;
1199  }
1200 
1201  mainN = fGdmlE->NewChild(nullptr, nullptr, "twistedtrap", nullptr);
1202  TString lname = GenName(geoShape->GetName(), TString::Format("%p", geoShape));
1203  fGdmlE->NewAttr(mainN, nullptr, "name", lname);
1204  if (IsNullParam(geoShape->GetDz(), "Dz", lname)) {
1205  return nullptr;
1206  }
1207 
1208  fGdmlE->NewAttr(mainN, nullptr, "z", TString::Format(fltPrecision.Data(), 2 * geoShape->GetDz()));
1209  fGdmlE->NewAttr(mainN, nullptr, "Theta", TString::Format(fltPrecision.Data(), geoShape->GetTheta()));
1210  fGdmlE->NewAttr(mainN, nullptr, "Phi", TString::Format(fltPrecision.Data(), geoShape->GetPhi()));
1211  fGdmlE->NewAttr(mainN, nullptr, "x1", TString::Format(fltPrecision.Data(), 2 * geoShape->GetBl1()));
1212  fGdmlE->NewAttr(mainN, nullptr, "x2", TString::Format(fltPrecision.Data(), 2 * geoShape->GetTl1()));
1213  fGdmlE->NewAttr(mainN, nullptr, "x3", TString::Format(fltPrecision.Data(), 2 * geoShape->GetBl2()));
1214  fGdmlE->NewAttr(mainN, nullptr, "x4", TString::Format(fltPrecision.Data(), 2 * geoShape->GetTl2()));
1215  fGdmlE->NewAttr(mainN, nullptr, "y1", TString::Format(fltPrecision.Data(), 2 * geoShape->GetH1()));
1216  fGdmlE->NewAttr(mainN, nullptr, "y2", TString::Format(fltPrecision.Data(), 2 * geoShape->GetH2()));
1217 
1218  fGdmlE->NewAttr(mainN, nullptr, "Alph", TString::Format(fltPrecision.Data(), geoShape->GetAlpha1()));
1219 
1220  //check if alpha1 equals to alpha2 (converting to string - to avoid problems with floats)
1221  if (TString::Format(fltPrecision.Data(), geoShape->GetAlpha1()) != TString::Format(fltPrecision.Data(), geoShape->GetAlpha2())) {
1222  Info("CreateTwistedTrapN",
1223  "ERROR! Object %s is not exported correctly because parameter Alpha2 is not declared in GDML schema",
1224  lname.Data());
1225  }
1226  //fGdmlE->NewAttr(mainN,0, "alpha2", TString::Format(fltPrecision.Data(), geoShape->GetAlpha2()));
1227  fGdmlE->NewAttr(mainN, nullptr, "PhiTwist", TString::Format(fltPrecision.Data(), geoShape->GetTwistAngle()));
1228 
1229  fGdmlE->NewAttr(mainN, nullptr, "aunit", "deg");
1230  fGdmlE->NewAttr(mainN, nullptr, "lunit", "cm");
1231  return mainN;
1232 }
1233 
1234 ////////////////////////////////////////////////////////////////////////////////
1235 /// Creates "trd" node for GDML from object TGeoTrd1
1236 
1237 XMLNodePointer_t TGDMLWrite::CreateTrdN(TGeoTrd1 * geoShape)
1238 {
1239  XMLNodePointer_t mainN = fGdmlE->NewChild(nullptr, nullptr, "trd", nullptr);
1240  const TString fltPrecision = TString::Format("%%.%dg", fFltPrecision);
1241  TString lname = GenName(geoShape->GetName(), TString::Format("%p", geoShape));
1242  fGdmlE->NewAttr(mainN, nullptr, "name", lname);
1243  if (IsNullParam(geoShape->GetDz(), "Dz", lname)) {
1244  return nullptr;
1245  }
1246 
1247  fGdmlE->NewAttr(mainN, nullptr, "x1", TString::Format(fltPrecision.Data(), 2 * geoShape->GetDx1()));
1248  fGdmlE->NewAttr(mainN, nullptr, "x2", TString::Format(fltPrecision.Data(), 2 * geoShape->GetDx2()));
1249  fGdmlE->NewAttr(mainN, nullptr, "y1", TString::Format(fltPrecision.Data(), 2 * geoShape->GetDy()));
1250  fGdmlE->NewAttr(mainN, nullptr, "y2", TString::Format(fltPrecision.Data(), 2 * geoShape->GetDy()));
1251  fGdmlE->NewAttr(mainN, nullptr, "z", TString::Format(fltPrecision.Data(), 2 * geoShape->GetDz()));
1252 
1253  fGdmlE->NewAttr(mainN, nullptr, "lunit", "cm");
1254  return mainN;
1255 }
1256 
1257 ////////////////////////////////////////////////////////////////////////////////
1258 /// Creates "trd" node for GDML from object TGeoTrd2
1259 
1260 XMLNodePointer_t TGDMLWrite::CreateTrdN(TGeoTrd2 * geoShape)
1261 {
1262  XMLNodePointer_t mainN = fGdmlE->NewChild(nullptr, nullptr, "trd", nullptr);
1263  const TString fltPrecision = TString::Format("%%.%dg", fFltPrecision);
1264  TString lname = GenName(geoShape->GetName(), TString::Format("%p", geoShape));
1265  fGdmlE->NewAttr(mainN, nullptr, "name", lname);
1266  if (IsNullParam(geoShape->GetDz(), "Dz", lname)) {
1267  return nullptr;
1268  }
1269 
1270  fGdmlE->NewAttr(mainN, nullptr, "x1", TString::Format(fltPrecision.Data(), 2 * geoShape->GetDx1()));
1271  fGdmlE->NewAttr(mainN, nullptr, "x2", TString::Format(fltPrecision.Data(), 2 * geoShape->GetDx2()));
1272  fGdmlE->NewAttr(mainN, nullptr, "y1", TString::Format(fltPrecision.Data(), 2 * geoShape->GetDy1()));
1273  fGdmlE->NewAttr(mainN, nullptr, "y2", TString::Format(fltPrecision.Data(), 2 * geoShape->GetDy2()));
1274  fGdmlE->NewAttr(mainN, nullptr, "z", TString::Format(fltPrecision.Data(), 2 * geoShape->GetDz()));
1275 
1276  fGdmlE->NewAttr(mainN, nullptr, "lunit", "cm");
1277  return mainN;
1278 }
1279 
1280 ////////////////////////////////////////////////////////////////////////////////
1281 /// Creates "tube" node for GDML from object TGeoTubeSeg
1282 
1283 XMLNodePointer_t TGDMLWrite::CreateTubeN(TGeoTubeSeg * geoShape)
1284 {
1285  XMLNodePointer_t mainN = fGdmlE->NewChild(nullptr, nullptr, "tube", nullptr);
1286  const TString fltPrecision = TString::Format("%%.%dg", fFltPrecision);
1287  TString lname = GenName(geoShape->GetName(), TString::Format("%p", geoShape));
1288  fGdmlE->NewAttr(mainN, nullptr, "name", lname);
1289  if (IsNullParam(geoShape->GetRmax(), "Rmax", lname) ||
1290  IsNullParam(geoShape->GetDz(), "Dz", lname)) {
1291  return nullptr;
1292  }
1293 
1294  fGdmlE->NewAttr(mainN, nullptr, "rmin", TString::Format(fltPrecision.Data(), geoShape->GetRmin()));
1295  fGdmlE->NewAttr(mainN, nullptr, "rmax", TString::Format(fltPrecision.Data(), geoShape->GetRmax()));
1296  fGdmlE->NewAttr(mainN, nullptr, "z", TString::Format(fltPrecision.Data(), 2 * geoShape->GetDz()));
1297  fGdmlE->NewAttr(mainN, nullptr, "startphi", TString::Format(fltPrecision.Data(), geoShape->GetPhi1()));
1298  fGdmlE->NewAttr(mainN, nullptr, "deltaphi", TString::Format(fltPrecision.Data(), geoShape->GetPhi2() - geoShape->GetPhi1()));
1299 
1300  fGdmlE->NewAttr(mainN, nullptr, "aunit", "deg");
1301  fGdmlE->NewAttr(mainN, nullptr, "lunit", "cm");
1302  return mainN;
1303 }
1304 
1305 ////////////////////////////////////////////////////////////////////////////////
1306 /// Creates "cutTube" node for GDML
1307 
1308 XMLNodePointer_t TGDMLWrite::CreateCutTubeN(TGeoCtub * geoShape)
1309 {
1310  XMLNodePointer_t mainN;
1311  const TString fltPrecision = TString::Format("%%.%dg", fFltPrecision);
1312 
1313  mainN = fGdmlE->NewChild(nullptr, nullptr, "cutTube", nullptr);
1314  TString lname = GenName(geoShape->GetName(), TString::Format("%p", geoShape));
1315  fGdmlE->NewAttr(mainN, nullptr, "name", lname);
1316  if (IsNullParam(geoShape->GetRmax(), "Rmax", lname) ||
1317  IsNullParam(geoShape->GetDz(), "Dz", lname)) {
1318  return nullptr;
1319  }
1320  //This is not needed, because cutTube is already supported by Geant4 9.5
1321  if (fgG4Compatibility == kTRUE && kFALSE) {
1322  TGeoShape * fakeCtub = CreateFakeCtub(geoShape);
1323  mainN = ChooseObject(fakeCtub);
1324 
1325  //register name for cuttube shape (so it will be found during volume export)
1326  lname = fNameList->fLst[TString::Format("%p", fakeCtub)];
1327  fNameList->fLst[TString::Format("%p", geoShape)] = lname;
1328  Info("CreateCutTubeN", "WARNING! %s - CutTube was replaced by intersection of TGeoTubSeg and two TGeoBBoxes",
1329  lname.Data());
1330  return mainN;
1331  }
1332  fGdmlE->NewAttr(mainN, nullptr, "rmin", TString::Format(fltPrecision.Data(), geoShape->GetRmin()));
1333  fGdmlE->NewAttr(mainN, nullptr, "rmax", TString::Format(fltPrecision.Data(), geoShape->GetRmax()));
1334  fGdmlE->NewAttr(mainN, nullptr, "z", TString::Format(fltPrecision.Data(), 2 * geoShape->GetDz()));
1335  fGdmlE->NewAttr(mainN, nullptr, "startphi", TString::Format(fltPrecision.Data(), geoShape->GetPhi1()));
1336  fGdmlE->NewAttr(mainN, nullptr, "deltaphi", TString::Format(fltPrecision.Data(), geoShape->GetPhi2() - geoShape->GetPhi1()));
1337  fGdmlE->NewAttr(mainN, nullptr, "lowX", TString::Format(fltPrecision.Data(), geoShape->GetNlow()[0]));
1338  fGdmlE->NewAttr(mainN, nullptr, "lowY", TString::Format(fltPrecision.Data(), geoShape->GetNlow()[1]));
1339  fGdmlE->NewAttr(mainN, nullptr, "lowZ", TString::Format(fltPrecision.Data(), geoShape->GetNlow()[2]));
1340  fGdmlE->NewAttr(mainN, nullptr, "highX", TString::Format(fltPrecision.Data(), geoShape->GetNhigh()[0]));
1341  fGdmlE->NewAttr(mainN, nullptr, "highY", TString::Format(fltPrecision.Data(), geoShape->GetNhigh()[1]));
1342  fGdmlE->NewAttr(mainN, nullptr, "highZ", TString::Format(fltPrecision.Data(), geoShape->GetNhigh()[2]));
1343 
1344  fGdmlE->NewAttr(mainN, nullptr, "aunit", "deg");
1345  fGdmlE->NewAttr(mainN, nullptr, "lunit", "cm");
1346 
1347  return mainN;
1348 }
1349 
1350 ////////////////////////////////////////////////////////////////////////////////
1351 /// Creates "tube" node for GDML from object TGeoTube
1352 
1353 XMLNodePointer_t TGDMLWrite::CreateTubeN(TGeoTube * geoShape)
1354 {
1355  XMLNodePointer_t mainN = fGdmlE->NewChild(nullptr, nullptr, "tube", nullptr);
1356  const TString fltPrecision = TString::Format("%%.%dg", fFltPrecision);
1357  TString lname = GenName(geoShape->GetName(), TString::Format("%p", geoShape));
1358  fGdmlE->NewAttr(mainN, nullptr, "name", lname);
1359  if (IsNullParam(geoShape->GetRmax(), "Rmax", lname) ||
1360  IsNullParam(geoShape->GetDz(), "Dz", lname)) {
1361  return nullptr;
1362  }
1363 
1364  fGdmlE->NewAttr(mainN, nullptr, "rmin", TString::Format(fltPrecision.Data(), geoShape->GetRmin()));
1365  fGdmlE->NewAttr(mainN, nullptr, "rmax", TString::Format(fltPrecision.Data(), geoShape->GetRmax()));
1366  fGdmlE->NewAttr(mainN, nullptr, "z", TString::Format(fltPrecision.Data(), 2 * geoShape->GetDz()));
1367  fGdmlE->NewAttr(mainN, nullptr, "startphi", TString::Format("%i", 0));
1368  fGdmlE->NewAttr(mainN, nullptr, "deltaphi", TString::Format("%i", 360));
1369 
1370  fGdmlE->NewAttr(mainN, nullptr, "aunit", "deg");
1371  fGdmlE->NewAttr(mainN, nullptr, "lunit", "cm");
1372  return mainN;
1373 }
1374 
1375 ////////////////////////////////////////////////////////////////////////////////
1376 /// Creates "zplane" node for GDML
1377 
1378 XMLNodePointer_t TGDMLWrite::CreateZplaneN(Double_t z, Double_t rmin, Double_t rmax)
1379 {
1380  XMLNodePointer_t mainN = fGdmlE->NewChild(nullptr, nullptr, "zplane", nullptr);
1381  const TString fltPrecision = TString::Format("%%.%dg", fFltPrecision);
1382 
1383  fGdmlE->NewAttr(mainN, nullptr, "z", TString::Format(fltPrecision.Data(), z));
1384  fGdmlE->NewAttr(mainN, nullptr, "rmin", TString::Format(fltPrecision.Data(), rmin));
1385  fGdmlE->NewAttr(mainN, nullptr, "rmax", TString::Format(fltPrecision.Data(), rmax));
1386 
1387  return mainN;
1388 }
1389 
1390 ////////////////////////////////////////////////////////////////////////////////
1391 /// Creates "polycone" node for GDML
1392 
1393 XMLNodePointer_t TGDMLWrite::CreatePolyconeN(TGeoPcon * geoShape)
1394 {
1395  XMLNodePointer_t mainN = fGdmlE->NewChild(nullptr, nullptr, "polycone", nullptr);
1396  const TString fltPrecision = TString::Format("%%.%dg", fFltPrecision);
1397  TString lname = GenName(geoShape->GetName(), TString::Format("%p", geoShape));
1398  fGdmlE->NewAttr(mainN, nullptr, "name", lname);
1399 
1400  fGdmlE->NewAttr(mainN, nullptr, "startphi", TString::Format(fltPrecision.Data(), geoShape->GetPhi1()));
1401  fGdmlE->NewAttr(mainN, nullptr, "deltaphi", TString::Format(fltPrecision.Data(), geoShape->GetDphi()));
1402 
1403  fGdmlE->NewAttr(mainN, nullptr, "aunit", "deg");
1404  fGdmlE->NewAttr(mainN, nullptr, "lunit", "cm");
1405  Int_t nZPlns = geoShape->GetNz();
1406  for (Int_t it = 0; it < nZPlns; it++) {
1407  //add zplane child node
1408  fGdmlE->AddChild(mainN, CreateZplaneN(geoShape->GetZ(it), geoShape->GetRmin(it), geoShape->GetRmax(it)));
1409  //compare actual plane and next plane
1410  if ((it < nZPlns - 1) && (geoShape->GetZ(it) == geoShape->GetZ(it + 1))) {
1411  //rmin of actual is greater then rmax of next one
1412  // | |rmax next
1413  // __ ...| |... __ < rmin actual
1414  // | | | |
1415  if (geoShape->GetRmin(it) > geoShape->GetRmax(it + 1)) {
1416  //adding plane from rmax next to rmin actual at the same z position
1417  if (fgG4Compatibility == kTRUE) {
1418  fGdmlE->AddChild(mainN, CreateZplaneN(geoShape->GetZ(it), geoShape->GetRmax(it + 1), geoShape->GetRmin(it)));
1419  Info("CreatePolyconeN", "WARNING! One plane was added to %s solid to be compatible with Geant4", lname.Data());
1420  } else {
1421  Info("CreatePolyconeN", "WARNING! Solid %s definition seems not contiguous may cause problems in Geant4", lname.Data());
1422  }
1423 
1424  }
1425  //rmin of next is greater then rmax of actual
1426  // | | | |
1427  // | |...___...| | rmin next
1428  // | | > rmax act
1429  if (geoShape->GetRmin(it + 1) > geoShape->GetRmax(it)) {
1430  //adding plane from rmax act to rmin next at the same z position
1431  if (fgG4Compatibility == kTRUE) {
1432  fGdmlE->AddChild(mainN, CreateZplaneN(geoShape->GetZ(it), geoShape->GetRmax(it), geoShape->GetRmin(it + 1)));
1433  Info("CreatePolyconeN", "WARNING! One plane was added to %s solid to be compatible with Geant4", lname.Data());
1434  } else {
1435  Info("CreatePolyconeN", "WARNING! Solid %s definition seems not contiguous may cause problems in Geant4", lname.Data());
1436  }
1437  }
1438  }
1439  }
1440  return mainN;
1441 }
1442 
1443 ////////////////////////////////////////////////////////////////////////////////
1444 /// Creates "torus" node for GDML
1445 
1446 XMLNodePointer_t TGDMLWrite::CreateTorusN(TGeoTorus * geoShape)
1447 {
1448  XMLNodePointer_t mainN = fGdmlE->NewChild(nullptr, nullptr, "torus", nullptr);
1449  const TString fltPrecision = TString::Format("%%.%dg", fFltPrecision);
1450  TString lname = GenName(geoShape->GetName(), TString::Format("%p", geoShape));
1451  fGdmlE->NewAttr(mainN, nullptr, "name", lname);
1452  if (IsNullParam(geoShape->GetRmax(), "Rmax", lname)) {
1453  return nullptr;
1454  }
1455 
1456  fGdmlE->NewAttr(mainN, nullptr, "rtor", TString::Format(fltPrecision.Data(), geoShape->GetR()));
1457  fGdmlE->NewAttr(mainN, nullptr, "rmin", TString::Format(fltPrecision.Data(), geoShape->GetRmin()));
1458  fGdmlE->NewAttr(mainN, nullptr, "rmax", TString::Format(fltPrecision.Data(), geoShape->GetRmax()));
1459  fGdmlE->NewAttr(mainN, nullptr, "startphi", TString::Format(fltPrecision.Data(), geoShape->GetPhi1()));
1460  fGdmlE->NewAttr(mainN, nullptr, "deltaphi", TString::Format(fltPrecision.Data(), geoShape->GetDphi()));
1461 
1462  fGdmlE->NewAttr(mainN, nullptr, "aunit", "deg");
1463  fGdmlE->NewAttr(mainN, nullptr, "lunit", "cm");
1464 
1465  return mainN;
1466 }
1467 
1468 ////////////////////////////////////////////////////////////////////////////////
1469 /// Creates "polyhedra" node for GDML
1470 
1471 XMLNodePointer_t TGDMLWrite::CreatePolyhedraN(TGeoPgon * geoShape)
1472 {
1473  XMLNodePointer_t mainN = fGdmlE->NewChild(nullptr, nullptr, "polyhedra", nullptr);
1474  const TString fltPrecision = TString::Format("%%.%dg", fFltPrecision);
1475  fGdmlE->NewAttr(mainN, nullptr, "name", GenName(geoShape->GetName(), TString::Format("%p", geoShape)));
1476 
1477  fGdmlE->NewAttr(mainN, nullptr, "startphi", TString::Format(fltPrecision.Data(), geoShape->GetPhi1()));
1478  fGdmlE->NewAttr(mainN, nullptr, "deltaphi", TString::Format(fltPrecision.Data(), geoShape->GetDphi()));
1479  fGdmlE->NewAttr(mainN, nullptr, "numsides", TString::Format("%i", geoShape->GetNedges()));
1480 
1481  fGdmlE->NewAttr(mainN, nullptr, "aunit", "deg");
1482  fGdmlE->NewAttr(mainN, nullptr, "lunit", "cm");
1483  for (Int_t it = 0; it < geoShape->GetNz(); it++) {
1484  //add zplane child node
1485  fGdmlE->AddChild(mainN, CreateZplaneN(geoShape->GetZ(it), geoShape->GetRmin(it), geoShape->GetRmax(it)));
1486  }
1487  return mainN;
1488 }
1489 
1490 ////////////////////////////////////////////////////////////////////////////////
1491 /// Creates "eltube" node for GDML
1492 
1493 XMLNodePointer_t TGDMLWrite::CreateEltubeN(TGeoEltu * geoShape)
1494 {
1495  XMLNodePointer_t mainN = fGdmlE->NewChild(nullptr, nullptr, "eltube", nullptr);
1496  const TString fltPrecision = TString::Format("%%.%dg", fFltPrecision);
1497  TString lname = GenName(geoShape->GetName(), TString::Format("%p", geoShape));
1498  fGdmlE->NewAttr(mainN, nullptr, "name", lname);
1499  if (IsNullParam(geoShape->GetA(), "A", lname) ||
1500  IsNullParam(geoShape->GetB(), "B", lname) ||
1501  IsNullParam(geoShape->GetDz(), "Dz", lname)) {
1502  return nullptr;
1503  }
1504 
1505  fGdmlE->NewAttr(mainN, nullptr, "dx", TString::Format(fltPrecision.Data(), geoShape->GetA()));
1506  fGdmlE->NewAttr(mainN, nullptr, "dy", TString::Format(fltPrecision.Data(), geoShape->GetB()));
1507  fGdmlE->NewAttr(mainN, nullptr, "dz", TString::Format(fltPrecision.Data(), geoShape->GetDz()));
1508 
1509  fGdmlE->NewAttr(mainN, nullptr, "lunit", "cm");
1510 
1511  return mainN;
1512 }
1513 
1514 ////////////////////////////////////////////////////////////////////////////////
1515 /// Creates "hype" node for GDML
1516 
1517 XMLNodePointer_t TGDMLWrite::CreateHypeN(TGeoHype * geoShape)
1518 {
1519  XMLNodePointer_t mainN = fGdmlE->NewChild(nullptr, nullptr, "hype", nullptr);
1520  const TString fltPrecision = TString::Format("%%.%dg", fFltPrecision);
1521  TString lname = GenName(geoShape->GetName(), TString::Format("%p", geoShape));
1522  fGdmlE->NewAttr(mainN, nullptr, "name", lname);
1523  if (IsNullParam(geoShape->GetDz(), "Dz", lname)) {
1524  return nullptr;
1525  }
1526 
1527 
1528  fGdmlE->NewAttr(mainN, nullptr, "rmin", TString::Format(fltPrecision.Data(), geoShape->GetRmin()));
1529  fGdmlE->NewAttr(mainN, nullptr, "rmax", TString::Format(fltPrecision.Data(), geoShape->GetRmax()));
1530  fGdmlE->NewAttr(mainN, nullptr, "inst", TString::Format(fltPrecision.Data(), geoShape->GetStIn()));
1531  fGdmlE->NewAttr(mainN, nullptr, "outst", TString::Format(fltPrecision.Data(), geoShape->GetStOut()));
1532  fGdmlE->NewAttr(mainN, nullptr, "z", TString::Format(fltPrecision.Data(), 2 * geoShape->GetDz()));
1533 
1534  fGdmlE->NewAttr(mainN, nullptr, "aunit", "deg");
1535  fGdmlE->NewAttr(mainN, nullptr, "lunit", "cm");
1536 
1537  return mainN;
1538 }
1539 
1540 ////////////////////////////////////////////////////////////////////////////////
1541 /// Creates "xtru" node for GDML
1542 
1543 XMLNodePointer_t TGDMLWrite::CreateXtrusionN(TGeoXtru * geoShape)
1544 {
1545  XMLNodePointer_t mainN = fGdmlE->NewChild(nullptr, nullptr, "xtru", nullptr);
1546  const TString fltPrecision = TString::Format("%%.%dg", fFltPrecision);
1547  TString lname = GenName(geoShape->GetName(), TString::Format("%p", geoShape));
1548  fGdmlE->NewAttr(mainN, nullptr, "name", lname);
1549 
1550  fGdmlE->NewAttr(mainN, nullptr, "lunit", "cm");
1551  XMLNodePointer_t childN;
1552  Int_t vertNum = geoShape->GetNvert();
1553  Int_t secNum = geoShape->GetNz();
1554  if (vertNum < 3 || secNum < 2) {
1555  Info("CreateXtrusionN", "ERROR! TGeoXtru %s has only %i vertices and %i sections. It was not exported",
1556  lname.Data(), vertNum, secNum);
1557  mainN = nullptr;
1558  return mainN;
1559  }
1560  for (Int_t it = 0; it < vertNum; it++) {
1561  //add twoDimVertex child node
1562  childN = fGdmlE->NewChild(nullptr, nullptr, "twoDimVertex", nullptr);
1563  fGdmlE->NewAttr(childN, nullptr, "x", TString::Format(fltPrecision.Data(), geoShape->GetX(it)));
1564  fGdmlE->NewAttr(childN, nullptr, "y", TString::Format(fltPrecision.Data(), geoShape->GetY(it)));
1565  fGdmlE->AddChild(mainN, childN);
1566  }
1567  for (Int_t it = 0; it < secNum; it++) {
1568  //add section child node
1569  childN = fGdmlE->NewChild(nullptr, nullptr, "section", nullptr);
1570  fGdmlE->NewAttr(childN, nullptr, "zOrder", TString::Format("%i", it));
1571  fGdmlE->NewAttr(childN, nullptr, "zPosition", TString::Format(fltPrecision.Data(), geoShape->GetZ(it)));
1572  fGdmlE->NewAttr(childN, nullptr, "xOffset", TString::Format(fltPrecision.Data(), geoShape->GetXOffset(it)));
1573  fGdmlE->NewAttr(childN, nullptr, "yOffset", TString::Format(fltPrecision.Data(), geoShape->GetYOffset(it)));
1574  fGdmlE->NewAttr(childN, nullptr, "scalingFactor", TString::Format(fltPrecision.Data(), geoShape->GetScale(it)));
1575  fGdmlE->AddChild(mainN, childN);
1576  }
1577  return mainN;
1578 }
1579 
1580 ////////////////////////////////////////////////////////////////////////////////
1581 /// Creates "ellipsoid" node for GDML
1582 /// this is a special case, because ellipsoid is not defined in ROOT
1583 /// so when intersection of scaled sphere and TGeoBBox is found,
1584 /// it is considered as an ellipsoid
1585 
1586 XMLNodePointer_t TGDMLWrite::CreateEllipsoidN(TGeoCompositeShape * geoShape, TString elName)
1587 {
1588  XMLNodePointer_t mainN = fGdmlE->NewChild(nullptr, nullptr, "ellipsoid", nullptr);
1589  const TString fltPrecision = TString::Format("%%.%dg", fFltPrecision);
1590  TGeoScaledShape *leftS = (TGeoScaledShape *)geoShape->GetBoolNode()->GetLeftShape(); //ScaledShape
1591  TGeoBBox *rightS = (TGeoBBox *)geoShape->GetBoolNode()->GetRightShape(); //BBox
1592 
1593 
1594  fGdmlE->NewAttr(mainN, nullptr, "name", elName.Data());
1595  Double_t sx = leftS->GetScale()->GetScale()[0];
1596  Double_t sy = leftS->GetScale()->GetScale()[1];
1597  Double_t radius = ((TGeoSphere *) leftS->GetShape())->GetRmax();
1598 
1599  Double_t ax, by, cz;
1600  cz = radius;
1601  ax = sx * radius;
1602  by = sy * radius;
1603 
1604  Double_t dz = rightS->GetDZ();
1605  Double_t zorig = rightS->GetOrigin()[2];
1606  Double_t zcut2 = dz + zorig;
1607  Double_t zcut1 = 2 * zorig - zcut2;
1608 
1609 
1610  fGdmlE->NewAttr(mainN, nullptr, "ax", TString::Format(fltPrecision.Data(), ax));
1611  fGdmlE->NewAttr(mainN, nullptr, "by", TString::Format(fltPrecision.Data(), by));
1612  fGdmlE->NewAttr(mainN, nullptr, "cz", TString::Format(fltPrecision.Data(), cz));
1613  fGdmlE->NewAttr(mainN, nullptr, "zcut1", TString::Format(fltPrecision.Data(), zcut1));
1614  fGdmlE->NewAttr(mainN, nullptr, "zcut2", TString::Format(fltPrecision.Data(), zcut2));
1615  fGdmlE->NewAttr(mainN, nullptr, "lunit", "cm");
1616 
1617  return mainN;
1618 }
1619 
1620 ////////////////////////////////////////////////////////////////////////////////
1621 /// Creates "elcone" (elliptical cone) node for GDML
1622 /// this is a special case, because elliptical cone is not defined in ROOT
1623 /// so when scaled cone is found, it is considered as a elliptical cone
1624 
1625 XMLNodePointer_t TGDMLWrite::CreateElConeN(TGeoScaledShape * geoShape)
1626 {
1627  XMLNodePointer_t mainN = fGdmlE->NewChild(nullptr, nullptr, "elcone", nullptr);
1628  const TString fltPrecision = TString::Format("%%.%dg", fFltPrecision);
1629  fGdmlE->NewAttr(mainN, nullptr, "name", GenName(geoShape->GetName(), TString::Format("%p", geoShape)));
1630  Double_t zcut = ((TGeoCone *) geoShape->GetShape())->GetDz();
1631  Double_t rx1 = ((TGeoCone *) geoShape->GetShape())->GetRmax1();
1632  Double_t rx2 = ((TGeoCone *) geoShape->GetShape())->GetRmax2();
1633  Double_t zmax = zcut * ((rx1 + rx2) / (rx1 - rx2));
1634  Double_t z = zcut + zmax;
1635 
1636  Double_t sy = geoShape->GetScale()->GetScale()[1];
1637  Double_t ry1 = sy * rx1;
1638 
1639  std::string format(TString::Format("%s/%s", fltPrecision.Data(), fltPrecision.Data()).Data());
1640  fGdmlE->NewAttr(mainN, nullptr, "dx", TString::Format(format.c_str(), rx1, z));
1641  fGdmlE->NewAttr(mainN, nullptr, "dy", TString::Format(format.c_str(), ry1, z));
1642  fGdmlE->NewAttr(mainN, nullptr, "zmax", TString::Format(fltPrecision.Data(), zmax));
1643  fGdmlE->NewAttr(mainN, nullptr, "zcut", TString::Format(fltPrecision.Data(), zcut));
1644  fGdmlE->NewAttr(mainN, nullptr, "lunit", "cm");
1645 
1646  return mainN;
1647 }
1648 
1649 ////////////////////////////////////////////////////////////////////////////////
1650 /// Creates common part of union intersection and subtraction nodes
1651 
1652 XMLNodePointer_t TGDMLWrite::CreateCommonBoolN(TGeoCompositeShape *geoShape)
1653 {
1654  XMLNodePointer_t mainN, ndR, ndL, childN;
1655  TString nodeName = GenName(geoShape->GetName(), TString::Format("%p", geoShape));
1656  TString lboolType;
1657  TGeoBoolNode::EGeoBoolType boolType = geoShape->GetBoolNode()->GetBooleanOperator();
1658  switch (boolType) {
1659  case TGeoBoolNode::kGeoUnion:
1660  lboolType = "union";
1661  break;
1662  case TGeoBoolNode::kGeoSubtraction:
1663  lboolType = "subtraction";
1664  break;
1665  case TGeoBoolNode::kGeoIntersection:
1666  lboolType = "intersection";
1667  break;
1668  }
1669 
1670  TGDMLWrite::Xyz lrot = GetXYZangles(geoShape->GetBoolNode()->GetLeftMatrix()->Inverse().GetRotationMatrix());
1671  const Double_t *ltr = geoShape->GetBoolNode()->GetLeftMatrix()->GetTranslation();
1672  TGDMLWrite::Xyz rrot = GetXYZangles(geoShape->GetBoolNode()->GetRightMatrix()->Inverse().GetRotationMatrix());
1673  const Double_t *rtr = geoShape->GetBoolNode()->GetRightMatrix()->GetTranslation();
1674 
1675  //specific case!
1676  //Ellipsoid tag preparing
1677  //if left == TGeoScaledShape AND right == TGeoBBox
1678  // AND if TGeoScaledShape->GetShape == TGeoSphere
1679  TGeoShape *leftS = geoShape->GetBoolNode()->GetLeftShape();
1680  TGeoShape *rightS = geoShape->GetBoolNode()->GetRightShape();
1681  if (strcmp(leftS->ClassName(), "TGeoScaledShape") == 0 &&
1682  strcmp(rightS->ClassName(), "TGeoBBox") == 0) {
1683  if (strcmp(((TGeoScaledShape *)leftS)->GetShape()->ClassName(), "TGeoSphere") == 0) {
1684  if (lboolType == "intersection") {
1685  mainN = CreateEllipsoidN(geoShape, nodeName);
1686  return mainN;
1687  }
1688  }
1689  }
1690 
1691  Xyz translL, translR;
1692  //translation
1693  translL.x = ltr[0];
1694  translL.y = ltr[1];
1695  translL.z = ltr[2];
1696  translR.x = rtr[0];
1697  translR.y = rtr[1];
1698  translR.z = rtr[2];
1699 
1700  //left and right nodes are created here also their names are created
1701  ndL = ChooseObject(geoShape->GetBoolNode()->GetLeftShape());
1702  ndR = ChooseObject(geoShape->GetBoolNode()->GetRightShape());
1703 
1704  //retrieve left and right node names by their pointer to make reference
1705  TString lname = fNameList->fLst[TString::Format("%p", geoShape->GetBoolNode()->GetLeftShape())];
1706  TString rname = fNameList->fLst[TString::Format("%p", geoShape->GetBoolNode()->GetRightShape())];
1707 
1708  //left and right nodes appended to main structure of nodes (if they are not already there)
1709  if (ndL != nullptr) {
1710  fGdmlE->AddChild(fSolidsNode, ndL);
1711  fSolCnt++;
1712  } else {
1713  if (lname.Contains("missing_") || lname == "") {
1714  Info("CreateCommonBoolN", "ERROR! Left node is NULL - Boolean Shape will be skipped");
1715  return nullptr;
1716  }
1717  }
1718  if (ndR != nullptr) {
1719  fGdmlE->AddChild(fSolidsNode, ndR);
1720  fSolCnt++;
1721  } else {
1722  if (rname.Contains("missing_") || rname == "") {
1723  Info("CreateCommonBoolN", "ERROR! Right node is NULL - Boolean Shape will be skipped");
1724  return nullptr;
1725  }
1726  }
1727 
1728  //create union node and its child nodes (or intersection or subtraction)
1729  /* <union name="...">
1730  * <first ref="left name" />
1731  * <second ref="right name" />
1732  * <firstposition .../>
1733  * <firstrotation .../>
1734  * <position .../>
1735  * <rotation .../>
1736  * </union>
1737  */
1738  mainN = fGdmlE->NewChild(nullptr, nullptr, lboolType.Data(), nullptr);
1739  fGdmlE->NewAttr(mainN, nullptr, "name", nodeName);
1740 
1741  //<first> (left)
1742  childN = fGdmlE->NewChild(nullptr, nullptr, "first", nullptr);
1743  fGdmlE->NewAttr(childN, nullptr, "ref", lname);
1744  fGdmlE->AddChild(mainN, childN);
1745 
1746  //<second> (right)
1747  childN = fGdmlE->NewChild(nullptr, nullptr, "second", nullptr);
1748  fGdmlE->NewAttr(childN, nullptr, "ref", rname);
1749  fGdmlE->AddChild(mainN, childN);
1750 
1751  //<firstposition> (left)
1752  if ((translL.x != 0.0) || (translL.y != 0.0) || (translL.z != 0.0)) {
1753  childN = CreatePositionN((nodeName + lname + "pos").Data(), translL, "firstposition");
1754  fGdmlE->AddChild(mainN, childN);
1755  }
1756  //<firstrotation> (left)
1757  if ((lrot.x != 0.0) || (lrot.y != 0.0) || (lrot.z != 0.0)) {
1758  childN = CreateRotationN((nodeName + lname + "rot").Data(), lrot, "firstrotation");
1759  fGdmlE->AddChild(mainN, childN);
1760  }
1761  //<position> (right)
1762  if ((translR.x != 0.0) || (translR.y != 0.0) || (translR.z != 0.0)) {
1763  childN = CreatePositionN((nodeName + rname + "pos").Data(), translR, "position");
1764  fGdmlE->AddChild(mainN, childN);
1765  }
1766  //<rotation> (right)
1767  if ((rrot.x != 0.0) || (rrot.y != 0.0) || (rrot.z != 0.0)) {
1768  childN = CreateRotationN((nodeName + rname + "rot").Data(), rrot, "rotation");
1769  fGdmlE->AddChild(mainN, childN);
1770  }
1771 
1772  return mainN;
1773 }
1774 
1775 ////////////////////////////////////////////////////////////////////////////////
1776 /// Creates "opticalsurface" node for GDML
1777 
1778 XMLNodePointer_t TGDMLWrite::CreateOpticalSurfaceN(TGeoOpticalSurface * geoSurf)
1779 {
1780  XMLNodePointer_t mainN = fGdmlE->NewChild(nullptr, nullptr, "opticalsurface", nullptr);
1781  const TString fltPrecision = TString::Format("%%.%dg", fFltPrecision);
1782  fGdmlE->NewAttr(mainN, nullptr, "name", geoSurf->GetName());
1783  fGdmlE->NewAttr(mainN, nullptr, "model", TGeoOpticalSurface::ModelToString(geoSurf->GetModel()));
1784  fGdmlE->NewAttr(mainN, nullptr, "finish", TGeoOpticalSurface::FinishToString(geoSurf->GetFinish()));
1785  fGdmlE->NewAttr(mainN, nullptr, "type", TGeoOpticalSurface::TypeToString(geoSurf->GetType()));
1786  fGdmlE->NewAttr(mainN, nullptr, "value", TString::Format(fltPrecision.Data(), geoSurf->GetValue()));
1787 
1788  // Write properties
1789  TList const &properties = geoSurf->GetProperties();
1790  if (properties.GetSize()) {
1791  TIter next(&properties);
1792  TNamed *property;
1793  while ((property = (TNamed*)next()))
1794  fGdmlE->AddChild(mainN, CreatePropertyN(*property));
1795  }
1796  return mainN;
1797 }
1798 
1799 ////////////////////////////////////////////////////////////////////////////////
1800 /// Creates "skinsurface" node for GDML
1801 
1802 XMLNodePointer_t TGDMLWrite::CreateSkinSurfaceN(TGeoSkinSurface * geoSurf)
1803 {
1804  XMLNodePointer_t mainN = fGdmlE->NewChild(nullptr, nullptr, "skinsurface", nullptr);
1805  fGdmlE->NewAttr(mainN, nullptr, "name", geoSurf->GetName());
1806  fGdmlE->NewAttr(mainN, nullptr, "surfaceproperty", geoSurf->GetTitle());
1807  // Cretate the logical volume reference node
1808  auto childN = fGdmlE->NewChild(nullptr, nullptr, "volumeref", nullptr);
1809  fGdmlE->NewAttr(childN, nullptr, "ref", geoSurf->GetVolume()->GetName());
1810  fGdmlE->AddChild(mainN, childN);
1811  return mainN;
1812 }
1813 
1814 ////////////////////////////////////////////////////////////////////////////////
1815 /// Creates "bordersurface" node for GDML
1816 
1817 XMLNodePointer_t TGDMLWrite::CreateBorderSurfaceN(TGeoBorderSurface * geoSurf)
1818 {
1819  XMLNodePointer_t mainN = fGdmlE->NewChild(nullptr, nullptr, "bordersurface", nullptr);
1820  fGdmlE->NewAttr(mainN, nullptr, "name", geoSurf->GetName());
1821  fGdmlE->NewAttr(mainN, nullptr, "surfaceproperty", geoSurf->GetTitle());
1822  // Cretate the logical volume reference node
1823  auto childN = fGdmlE->NewChild(nullptr, nullptr, "physvolref", nullptr);
1824  fGdmlE->NewAttr(childN, nullptr, "ref", geoSurf->GetNode1()->GetName());
1825  fGdmlE->NewAttr(childN, nullptr, "ref", geoSurf->GetNode2()->GetName());
1826  fGdmlE->AddChild(mainN, childN);
1827  return mainN;
1828 }
1829 
1830 ////////////////////////////////////////////////////////////////////////////////
1831 /// Creates "position" kind of node for GDML
1832 
1833 XMLNodePointer_t TGDMLWrite::CreatePositionN(const char * name, Xyz position, const char * type, const char * unit)
1834 {
1835  XMLNodePointer_t mainN = fGdmlE->NewChild(nullptr, nullptr, type, nullptr);
1836  const TString fltPrecision = TString::Format("%%.%dg", fFltPrecision);
1837  fGdmlE->NewAttr(mainN, nullptr, "name", name);
1838  fGdmlE->NewAttr(mainN, nullptr, "x", TString::Format(fltPrecision.Data(), position.x));
1839  fGdmlE->NewAttr(mainN, nullptr, "y", TString::Format(fltPrecision.Data(), position.y));
1840  fGdmlE->NewAttr(mainN, nullptr, "z", TString::Format(fltPrecision.Data(), position.z));
1841  fGdmlE->NewAttr(mainN, nullptr, "unit", unit);
1842  return mainN;
1843 }
1844 
1845 ////////////////////////////////////////////////////////////////////////////////
1846 /// Creates "rotation" kind of node for GDML
1847 
1848 XMLNodePointer_t TGDMLWrite::CreateRotationN(const char * name, Xyz rotation, const char * type, const char * unit)
1849 {
1850  XMLNodePointer_t mainN = fGdmlE->NewChild(nullptr, nullptr, type, nullptr);
1851  const TString fltPrecision = TString::Format("%%.%dg", fFltPrecision);
1852  fGdmlE->NewAttr(mainN, nullptr, "name", name);
1853  fGdmlE->NewAttr(mainN, nullptr, "x", TString::Format(fltPrecision.Data(), rotation.x));
1854  fGdmlE->NewAttr(mainN, nullptr, "y", TString::Format(fltPrecision.Data(), rotation.y));
1855  fGdmlE->NewAttr(mainN, nullptr, "z", TString::Format(fltPrecision.Data(), rotation.z));
1856  fGdmlE->NewAttr(mainN, nullptr, "unit", unit);
1857  return mainN;
1858 }
1859 
1860 ////////////////////////////////////////////////////////////////////////////////
1861 /// Creates "matrix" kind of node for GDML
1862 
1863 XMLNodePointer_t TGDMLWrite::CreateMatrixN(TGDMLMatrix const *matrix)
1864 {
1865  std::stringstream vals;
1866  size_t cols = matrix->GetCols();
1867  size_t rows = matrix->GetRows();
1868  XMLNodePointer_t mainN = fGdmlE->NewChild(nullptr, nullptr, "matrix", nullptr);
1869  fGdmlE->NewAttr(mainN, nullptr, "name", matrix->GetName());
1870  fGdmlE->NewAttr(mainN, nullptr, "coldim", TString::Format("%zu", cols));
1871  for(size_t i=0; i<rows; ++i) {
1872  for(size_t j=0; j<cols; ++j) {
1873  vals << matrix->Get(i,j);
1874  if ( j < cols-1 ) vals << ' ';
1875  }
1876  if ( i < rows-1 ) vals << '\n';
1877  }
1878  fGdmlE->NewAttr(mainN, nullptr, "values", vals.str().c_str());
1879  return mainN;
1880 }
1881 
1882 ////////////////////////////////////////////////////////////////////////////////
1883 /// Creates "constant" kind of node for GDML
1884 
1885 XMLNodePointer_t TGDMLWrite::CreateConstantN(const char *name, Double_t value)
1886 {
1887  XMLNodePointer_t mainN = fGdmlE->NewChild(nullptr, nullptr, "constant", nullptr);
1888  const TString fltPrecision = TString::Format("%%.%dg", fFltPrecision);
1889  fGdmlE->NewAttr(mainN, nullptr, "name", name);
1890  fGdmlE->NewAttr(mainN, nullptr, "value", TString::Format(fltPrecision.Data(), value));
1891  return mainN;
1892 }
1893 
1894 ////////////////////////////////////////////////////////////////////////////////
1895 /// Creates "setup" node for GDML
1896 
1897 XMLNodePointer_t TGDMLWrite::CreateSetupN(const char * topVolName, const char * name, const char * version)
1898 {
1899  XMLNodePointer_t setupN = fGdmlE->NewChild(nullptr, nullptr, "setup", nullptr);
1900  fGdmlE->NewAttr(setupN, nullptr, "name", name);
1901  fGdmlE->NewAttr(setupN, nullptr, "version", version);
1902  XMLNodePointer_t fworldN = fGdmlE->NewChild(setupN, nullptr, "world", nullptr);
1903  fGdmlE->NewAttr(fworldN, nullptr, "ref", topVolName);
1904  return setupN;
1905 }
1906 
1907 ////////////////////////////////////////////////////////////////////////////////
1908 /// Creates "volume" node for GDML
1909 
1910 XMLNodePointer_t TGDMLWrite::StartVolumeN(const char * name, const char * solid, const char * material)
1911 {
1912  XMLNodePointer_t childN;
1913  XMLNodePointer_t mainN = fGdmlE->NewChild(nullptr, nullptr, "volume", nullptr);
1914  fGdmlE->NewAttr(mainN, nullptr, "name", name);
1915 
1916  childN = fGdmlE->NewChild(nullptr, nullptr, "materialref", nullptr);
1917  fGdmlE->NewAttr(childN, nullptr, "ref", material);
1918  fGdmlE->AddChild(mainN, childN);
1919 
1920  childN = fGdmlE->NewChild(nullptr, nullptr, "solidref", nullptr);
1921  fGdmlE->NewAttr(childN, nullptr, "ref", solid);
1922  fGdmlE->AddChild(mainN, childN);
1923 
1924  return mainN;
1925 }
1926 
1927 ////////////////////////////////////////////////////////////////////////////////
1928 /// Creates "assembly" node for GDML
1929 
1930 XMLNodePointer_t TGDMLWrite::StartAssemblyN(const char * name)
1931 {
1932  XMLNodePointer_t mainN = fGdmlE->NewChild(nullptr, nullptr, "assembly", nullptr);
1933  fGdmlE->NewAttr(mainN, nullptr, "name", name);
1934 
1935  return mainN;
1936 }
1937 
1938 ////////////////////////////////////////////////////////////////////////////////
1939 /// Creates "physvol" node for GDML
1940 
1941 XMLNodePointer_t TGDMLWrite::CreatePhysVolN(const char *name, Int_t copyno, const char * volref, const char * posref, const char * rotref, XMLNodePointer_t scaleN)
1942 {
1943  fPhysVolCnt++;
1944  XMLNodePointer_t childN;
1945  XMLNodePointer_t mainN = fGdmlE->NewChild(nullptr, nullptr, "physvol", nullptr);
1946  fGdmlE->NewAttr(mainN, nullptr, "name", name);
1947  fGdmlE->NewAttr(mainN, nullptr, "copynumber", TString::Format("%d",copyno));
1948 
1949  childN = fGdmlE->NewChild(nullptr, nullptr, "volumeref", nullptr);
1950  fGdmlE->NewAttr(childN, nullptr, "ref", volref);
1951  fGdmlE->AddChild(mainN, childN);
1952 
1953  childN = fGdmlE->NewChild(nullptr, nullptr, "positionref", nullptr);
1954  fGdmlE->NewAttr(childN, nullptr, "ref", posref);
1955  fGdmlE->AddChild(mainN, childN);
1956 
1957  //if is not empty string add this node
1958  if (strcmp(rotref, "") != 0) {
1959  childN = fGdmlE->NewChild(nullptr, nullptr, "rotationref", nullptr);
1960  fGdmlE->NewAttr(childN, nullptr, "ref", rotref);
1961  fGdmlE->AddChild(mainN, childN);
1962  }
1963  if (scaleN != nullptr) {
1964  fGdmlE->AddChild(mainN, scaleN);
1965  }
1966 
1967  return mainN;
1968 }
1969 
1970 ////////////////////////////////////////////////////////////////////////////////
1971 /// Creates "divisionvol" node for GDML
1972 
1973 XMLNodePointer_t TGDMLWrite::CreateDivisionN(Double_t offset, Double_t width, Int_t number, const char * axis, const char * unit, const char * volref)
1974 {
1975  XMLNodePointer_t childN = 0;
1976  XMLNodePointer_t mainN = fGdmlE->NewChild(nullptr, nullptr, "divisionvol", nullptr);
1977  fGdmlE->NewAttr(mainN, nullptr, "axis", axis);
1978  fGdmlE->NewAttr(mainN, nullptr, "number", TString::Format("%i", number));
1979  const TString fltPrecision = TString::Format("%%.%dg", fFltPrecision);
1980  if (fgG4Compatibility == kTRUE) {
1981  //if eg. full length is 20 and width * number = 20,0001 problem in geant4
1982  //unit is either in cm or degrees nothing else
1983  width = (floor(width * 1E4)) * 1E-4;
1984  if ((offset >= 0.) && (strcmp(axis, "kPhi") == 0)) {
1985  Int_t offsetI = (Int_t) offset;
1986  Double_t decimals = offset - offsetI;
1987  //put to range from 0 to 360 add decimals and then put to range 0 -> -360
1988  offset = (offsetI % 360) + decimals - 360;
1989  }
1990  }
1991  fGdmlE->NewAttr(mainN, nullptr, "width", TString::Format(fltPrecision.Data(), width));
1992 
1993  fGdmlE->NewAttr(mainN, nullptr, "offset", TString::Format(fltPrecision.Data(), offset));
1994  fGdmlE->NewAttr(mainN, nullptr, "unit", unit);
1995  if (strcmp(volref, "") != 0) {
1996  childN = fGdmlE->NewChild(nullptr, nullptr, "volumeref", nullptr);
1997  fGdmlE->NewAttr(childN, nullptr, "ref", volref);
1998  }
1999  fGdmlE->AddChild(mainN, childN);
2000 
2001 
2002  return mainN;
2003 }
2004 
2005 ////////////////////////////////////////////////////////////////////////////////
2006 /// Chooses the object and method that should be used for processing object
2007 
2008 XMLNodePointer_t TGDMLWrite::ChooseObject(TGeoShape *geoShape)
2009 {
2010  const char * clsname = geoShape->ClassName();
2011  XMLNodePointer_t solidN;
2012 
2013  if (CanProcess((TObject *)geoShape) == kFALSE) {
2014  return nullptr;
2015  }
2016 
2017  //process different shapes
2018  if (strcmp(clsname, "TGeoBBox") == 0) {
2019  solidN = CreateBoxN((TGeoBBox*) geoShape);
2020  } else if (strcmp(clsname, "TGeoParaboloid") == 0) {
2021  solidN = CreateParaboloidN((TGeoParaboloid*) geoShape);
2022  } else if (strcmp(clsname, "TGeoSphere") == 0) {
2023  solidN = CreateSphereN((TGeoSphere*) geoShape);
2024  } else if (strcmp(clsname, "TGeoArb8") == 0) {
2025  solidN = CreateArb8N((TGeoArb8*) geoShape);
2026  } else if (strcmp(clsname, "TGeoConeSeg") == 0) {
2027  solidN = CreateConeN((TGeoConeSeg*) geoShape);
2028  } else if (strcmp(clsname, "TGeoCone") == 0) {
2029  solidN = CreateConeN((TGeoCone*) geoShape);
2030  } else if (strcmp(clsname, "TGeoPara") == 0) {
2031  solidN = CreateParaN((TGeoPara*) geoShape);
2032  } else if (strcmp(clsname, "TGeoTrap") == 0) {
2033  solidN = CreateTrapN((TGeoTrap*) geoShape);
2034  } else if (strcmp(clsname, "TGeoGtra") == 0) {
2035  solidN = CreateTwistedTrapN((TGeoGtra*) geoShape);
2036  } else if (strcmp(clsname, "TGeoTrd1") == 0) {
2037  solidN = CreateTrdN((TGeoTrd1*) geoShape);
2038  } else if (strcmp(clsname, "TGeoTrd2") == 0) {
2039  solidN = CreateTrdN((TGeoTrd2*) geoShape);
2040  } else if (strcmp(clsname, "TGeoTubeSeg") == 0) {
2041  solidN = CreateTubeN((TGeoTubeSeg*) geoShape);
2042  } else if (strcmp(clsname, "TGeoCtub") == 0) {
2043  solidN = CreateCutTubeN((TGeoCtub*) geoShape);
2044  } else if (strcmp(clsname, "TGeoTube") == 0) {
2045  solidN = CreateTubeN((TGeoTube*) geoShape);
2046  } else if (strcmp(clsname, "TGeoPcon") == 0) {
2047  solidN = CreatePolyconeN((TGeoPcon*) geoShape);
2048  } else if (strcmp(clsname, "TGeoTorus") == 0) {
2049  solidN = CreateTorusN((TGeoTorus*) geoShape);
2050  } else if (strcmp(clsname, "TGeoPgon") == 0) {
2051  solidN = CreatePolyhedraN((TGeoPgon*) geoShape);
2052  } else if (strcmp(clsname, "TGeoEltu") == 0) {
2053  solidN = CreateEltubeN((TGeoEltu*) geoShape);
2054  } else if (strcmp(clsname, "TGeoHype") == 0) {
2055  solidN = CreateHypeN((TGeoHype*) geoShape);
2056  } else if (strcmp(clsname, "TGeoXtru") == 0) {
2057  solidN = CreateXtrusionN((TGeoXtru*) geoShape);
2058  } else if (strcmp(clsname, "TGeoScaledShape") == 0) {
2059  TGeoScaledShape * geoscale = (TGeoScaledShape *) geoShape;
2060  TString scaleObjClsName = geoscale->GetShape()->ClassName();
2061  if (scaleObjClsName == "TGeoCone") {
2062  solidN = CreateElConeN((TGeoScaledShape*) geoShape);
2063  } else {
2064  Info("ChooseObject",
2065  "ERROR! TGeoScaledShape object is not possible to process correctly. %s object is processed without scale",
2066  scaleObjClsName.Data());
2067  solidN = ChooseObject(geoscale->GetShape());
2068  //Name has to be propagated to geoscale level pointer
2069  fNameList->fLst[TString::Format("%p", geoscale)] =
2070  fNameList->fLst[TString::Format("%p", geoscale->GetShape())];
2071  }
2072  } else if (strcmp(clsname, "TGeoCompositeShape") == 0) {
2073  solidN = CreateCommonBoolN((TGeoCompositeShape*) geoShape);
2074  } else if (strcmp(clsname, "TGeoUnion") == 0) {
2075  solidN = CreateCommonBoolN((TGeoCompositeShape*) geoShape);
2076  } else if (strcmp(clsname, "TGeoIntersection") == 0) {
2077  solidN = CreateCommonBoolN((TGeoCompositeShape*) geoShape);
2078  } else if (strcmp(clsname, "TGeoSubtraction") == 0) {
2079  solidN = CreateCommonBoolN((TGeoCompositeShape*) geoShape);
2080  } else {
2081  Info("ChooseObject", "ERROR! %s Solid CANNOT be processed, solid is NOT supported",
2082  clsname);
2083  solidN = nullptr;
2084  }
2085  if (solidN == nullptr) {
2086  if (fNameList->fLst[TString::Format("%p", geoShape)] == "") {
2087  TString missingName = geoShape->GetName();
2088  GenName("missing_" + missingName, TString::Format("%p", geoShape));
2089  } else {
2090  fNameList->fLst[TString::Format("%p", geoShape)] = "missing_" + fNameList->fLst[TString::Format("%p", geoShape)];
2091  }
2092  }
2093 
2094  return solidN;
2095 }
2096 
2097 ////////////////////////////////////////////////////////////////////////////////
2098 /// Retrieves X Y Z angles from rotation matrix
2099 
2100 TGDMLWrite::Xyz TGDMLWrite::GetXYZangles(const Double_t * rotationMatrix)
2101 {
2102  TGDMLWrite::Xyz lxyz;
2103  Double_t a, b, c;
2104  Double_t rad = 180.0 / TMath::ACos(-1.0);
2105  const Double_t *r = rotationMatrix;
2106  Double_t cosb = TMath::Sqrt(r[0] * r[0] + r[1] * r[1]);
2107  if (cosb > 0.00001) {
2108  a = TMath::ATan2(r[5], r[8]) * rad;
2109  b = TMath::ATan2(-r[2], cosb) * rad;
2110  c = TMath::ATan2(r[1], r[0]) * rad;
2111  } else {
2112  a = TMath::ATan2(-r[7], r[4]) * rad;
2113  b = TMath::ATan2(-r[2], cosb) * rad;
2114  c = 0;
2115  }
2116  lxyz.x = a;
2117  lxyz.y = b;
2118  lxyz.z = c;
2119  return lxyz;
2120 }
2121 
2122 ////////////////////////////////////////////////////////////////////////////////
2123 /// Method creating cutTube as an intersection of tube and two boxes
2124 /// - not used anymore because cutube is supported in Geant4 9.5
2125 
2126 TGeoCompositeShape* TGDMLWrite::CreateFakeCtub(TGeoCtub* geoShape)
2127 {
2128  Double_t rmin = geoShape->GetRmin();
2129  Double_t rmax = geoShape->GetRmax();
2130  Double_t z = geoShape->GetDz();
2131  Double_t startphi = geoShape->GetPhi1();
2132  Double_t deltaphi = geoShape->GetPhi2();
2133  Double_t x1 = geoShape->GetNlow()[0];
2134  Double_t y1 = geoShape->GetNlow()[1];
2135  Double_t z1 = geoShape->GetNlow()[2];
2136  Double_t x2 = geoShape->GetNhigh()[0];
2137  Double_t y2 = geoShape->GetNhigh()[1];
2138  Double_t z2 = geoShape->GetNhigh()[2];
2139  TString xname = geoShape->GetName();
2140 
2141 
2142  Double_t h0 = 2.*((TGeoBBox*)geoShape)->GetDZ();
2143  Double_t h1 = 2 * z;
2144  Double_t h2 = 2 * z;
2145  Double_t boxdx = 1E8 * (2 * rmax) + (2 * z);
2146 
2147  TGeoTubeSeg *T = new TGeoTubeSeg((xname + "T").Data(), rmin, rmax, h0, startphi, deltaphi);
2148  TGeoBBox *B1 = new TGeoBBox((xname + "B1").Data(), boxdx, boxdx, h1);
2149  TGeoBBox *B2 = new TGeoBBox((xname + "B2").Data(), boxdx, boxdx, h2);
2150 
2151 
2152  //first box position parameters
2153  Double_t phi1 = 360 - TMath::ATan2(x1, y1) * TMath::RadToDeg();
2154  Double_t theta1 = 360 - TMath::ATan2(sqrt(x1 * x1 + y1 * y1), z1) * TMath::RadToDeg();
2155 
2156  Double_t phi11 = TMath::ATan2(y1, x1) * TMath::RadToDeg() ;
2157  Double_t theta11 = TMath::ATan2(z1, sqrt(x1 * x1 + y1 * y1)) * TMath::RadToDeg() ;
2158 
2159  Double_t xpos1 = h1 * TMath::Cos((theta11) * TMath::DegToRad()) * TMath::Cos((phi11) * TMath::DegToRad()) * (-1);
2160  Double_t ypos1 = h1 * TMath::Cos((theta11) * TMath::DegToRad()) * TMath::Sin((phi11) * TMath::DegToRad()) * (-1);
2161  Double_t zpos1 = h1 * TMath::Sin((theta11) * TMath::DegToRad()) * (-1);
2162 
2163  //second box position parameters
2164  Double_t phi2 = 360 - TMath::ATan2(x2, y2) * TMath::RadToDeg();
2165  Double_t theta2 = 360 - TMath::ATan2(sqrt(x2 * x2 + y2 * y2), z2) * TMath::RadToDeg();
2166 
2167  Double_t phi21 = TMath::ATan2(y2, x2) * TMath::RadToDeg() ;
2168  Double_t theta21 = TMath::ATan2(z2, sqrt(x2 * x2 + y2 * y2)) * TMath::RadToDeg() ;
2169 
2170  Double_t xpos2 = h2 * TMath::Cos((theta21) * TMath::DegToRad()) * TMath::Cos((phi21) * TMath::DegToRad()) * (-1);
2171  Double_t ypos2 = h2 * TMath::Cos((theta21) * TMath::DegToRad()) * TMath::Sin((phi21) * TMath::DegToRad()) * (-1);
2172  Double_t zpos2 = h2 * TMath::Sin((theta21) * TMath::DegToRad()) * (-1);
2173 
2174 
2175  //positioning
2176  TGeoTranslation *t0 = new TGeoTranslation(0, 0, 0);
2177  TGeoTranslation *t1 = new TGeoTranslation(0 + xpos1, 0 + ypos1, 0 + (zpos1 - z));
2178  TGeoTranslation *t2 = new TGeoTranslation(0 + xpos2, 0 + ypos2, 0 + (zpos2 + z));
2179  TGeoRotation *r0 = new TGeoRotation((xname + "r0").Data());
2180  TGeoRotation *r1 = new TGeoRotation((xname + "r1").Data());
2181  TGeoRotation *r2 = new TGeoRotation((xname + "r2").Data());
2182 
2183  r1->SetAngles(phi1, theta1, 0);
2184  r2->SetAngles(phi2, theta2, 0);
2185 
2186  TGeoMatrix* m0 = new TGeoCombiTrans(*t0, *r0);
2187  TGeoMatrix* m1 = new TGeoCombiTrans(*t1, *r1);
2188  TGeoMatrix* m2 = new TGeoCombiTrans(*t2, *r2);
2189 
2190  TGeoCompositeShape *CS1 = new TGeoCompositeShape((xname + "CS1").Data(), new TGeoIntersection(T, B1, m0, m1));
2191  TGeoCompositeShape *cs = new TGeoCompositeShape((xname + "CS").Data(), new TGeoIntersection(CS1, B2, m0, m2));
2192  delete t0;
2193  delete t1;
2194  delete t2;
2195  delete r0;
2196  delete r1;
2197  delete r2;
2198  return cs;
2199 }
2200 
2201 ////////////////////////////////////////////////////////////////////////////////
2202 /// Checks whether name2check is in (NameList) list
2203 
2204 Bool_t TGDMLWrite::IsInList(NameList list, TString name2check)
2205 {
2206  Bool_t isIN = list[name2check];
2207  return isIN;
2208 }
2209 
2210 ////////////////////////////////////////////////////////////////////////////////
2211 ///NCNAME basic restrictions
2212 ///Replace "$" character with empty character etc.
2213 
2214 TString TGDMLWrite::GenName(TString oldname)
2215 {
2216  TString newname = oldname.ReplaceAll("$", "");
2217  newname = newname.ReplaceAll(" ", "_");
2218  // :, @, $, %, &, /, +, ,, ;, whitespace characters or different parenthesis
2219  newname = newname.ReplaceAll(":", "");
2220  newname = newname.ReplaceAll("@", "");
2221  newname = newname.ReplaceAll("%", "");
2222  newname = newname.ReplaceAll("&", "");
2223  newname = newname.ReplaceAll("/", "");
2224  newname = newname.ReplaceAll("+", "");
2225  newname = newname.ReplaceAll(";", "");
2226  newname = newname.ReplaceAll("{", "");
2227  newname = newname.ReplaceAll("}", "");
2228  newname = newname.ReplaceAll("(", "");
2229  newname = newname.ReplaceAll(")", "");
2230  newname = newname.ReplaceAll("[", "");
2231  newname = newname.ReplaceAll("]", "");
2232  newname = newname.ReplaceAll("_refl", "");
2233  //workaround if first letter is digit than replace it to "O" (ou character)
2234  TString fstLet = newname(0, 1);
2235  if (fstLet.IsDigit()) {
2236  newname = "O" + newname(1, newname.Length());
2237  }
2238  return newname;
2239 }
2240 
2241 ////////////////////////////////////////////////////////////////////////////////
2242 /// Important function which is responsible for naming volumes, solids and materials
2243 
2244 TString TGDMLWrite::GenName(TString oldname, TString objPointer)
2245 {
2246  TString newname = GenName(oldname);
2247  if (newname != oldname) {
2248  if (fgkMaxNameErr > fActNameErr) {
2249  Info("GenName",
2250  "WARNING! Name of the object was changed because it failed to comply with NCNAME xml datatype restrictions.");
2251  } else if ((fgkMaxNameErr == fActNameErr)) {
2252  Info("GenName",
2253  "WARNING! Probably more names are going to be changed to comply with NCNAME xml datatype restriction, but it will not be displayed on the screen.");
2254  }
2255  fActNameErr++;
2256  }
2257  TString nameIter;
2258  Int_t iter = 0;
2259  switch (fgNamingSpeed) {
2260  case kfastButUglySufix:
2261  newname = newname + "0x" + objPointer;
2262  break;
2263  case kelegantButSlow:
2264  //0 means not in the list
2265  iter = fNameList->fLstIter[newname];
2266  if (iter == 0) {
2267  nameIter = "";
2268  } else {
2269  nameIter = TString::Format("0x%i", iter);
2270  }
2271  fNameList->fLstIter[newname]++;
2272  newname = newname + nameIter;
2273  break;
2274  case kwithoutSufixNotUniq:
2275  //no change
2276  break;
2277  }
2278  //store the name (mapped to pointer)
2279  fNameList->fLst[objPointer] = newname;
2280  return newname;
2281 }
2282 
2283 
2284 ////////////////////////////////////////////////////////////////////////////////
2285 /// Method which tests whether solids can be processed
2286 
2287 Bool_t TGDMLWrite::CanProcess(TObject *pointer)
2288 {
2289  Bool_t isProcessed = kFALSE;
2290  isProcessed = pointer->TestBit(fgkProcBit);
2291  pointer->SetBit(fgkProcBit, kTRUE);
2292  return !(isProcessed);
2293 }
2294 
2295 ////////////////////////////////////////////////////////////////////////////////
2296 /// Method that retrieves axis and unit along which object is divided
2297 
2298 TString TGDMLWrite::GetPattAxis(Int_t divAxis, const char * pattName, TString& unit)
2299 {
2300  TString resaxis;
2301  unit = "cm";
2302  switch (divAxis) {
2303  case 1:
2304  if (strcmp(pattName, "TGeoPatternX") == 0) {
2305  return "kXAxis";
2306  } else if (strcmp(pattName, "TGeoPatternCylR") == 0) {
2307  return "kRho";
2308  }
2309  break;
2310  case 2:
2311  if (strcmp(pattName, "TGeoPatternY") == 0) {
2312  return "kYAxis";
2313  } else if (strcmp(pattName, "TGeoPatternCylPhi") == 0) {
2314  unit = "deg";
2315  return "kPhi";
2316  }
2317  break;
2318  case 3:
2319  if (strcmp(pattName, "TGeoPatternZ") == 0) {
2320  return "kZAxis";
2321  }
2322  break;
2323  default:
2324  return "kUndefined";
2325  break;
2326  }
2327  return "kUndefined";
2328 }
2329 
2330 ////////////////////////////////////////////////////////////////////////////////
2331 /// Check for null parameter to skip the NULL objects
2332 
2333 Bool_t TGDMLWrite::IsNullParam(Double_t parValue, TString parName, TString objName)
2334 {
2335  if (parValue == 0.) {
2336  Info("IsNullParam", "ERROR! %s is NULL due to %s = %.12g, Volume based on this shape will be skipped",
2337  objName.Data(),
2338  parName.Data(),
2339  parValue);
2340  return kTRUE;
2341  }
2342  return kFALSE;
2343 }
2344 
2345 ////////////////////////////////////////////////////////////////////////////////
2346 /// Unsetting bits that were changed in gGeoManager during export so that export
2347 /// can be run more times with the same instance of gGeoManager.
2348 
2349 void TGDMLWrite::UnsetTemporaryBits(TGeoManager * geoMng)
2350 {
2351  TIter next(geoMng->GetListOfVolumes());
2352  TGeoVolume *vol;
2353  while ((vol = (TGeoVolume *)next())) {
2354  ((TObject *)vol->GetShape())->SetBit(fgkProcBit, kFALSE);
2355  vol->SetAttBit(fgkProcBitVol, kFALSE);
2356  }
2357 
2358 }
2359 
2360 
2361 ////////////////////////////////////////////////////////////////////////////////
2362 //
2363 // Backwards compatibility for old DD4hep version (to be removed in the future)
2364 //
2365 ////////////////////////////////////////////////////////////////////////////////
2366 
2367 ////////////////////////////////////////////////////////////////////////////////
2368 // Backwards compatibility (to be removed in the future): Wrapper to only selectively write one branch
2369 void TGDMLWrite::WriteGDMLfile(TGeoManager * geomanager, TGeoVolume* volume, const char* filename, TString option)
2370 {
2371  TList materials, volumes, nodes;
2372  MaterialExtractor extract;
2373  if ( !volume ) {
2374  Info("WriteGDMLfile", "Invalid Volume reference to extract GDML information!");
2375  return;
2376  }
2377  extract(volume);
2378  for(TGeoMaterial* m : extract.materials)
2379  materials.Add(m);
2380  fTopVolumeName = volume->GetName();
2381  fSurfaceList.clear();
2382  fVolumeList.clear();
2383  fNodeList.clear();
2384  WriteGDMLfile(geomanager, volume, &materials, filename, option);
2385  materials.Clear("nodelete");
2386  volumes.Clear("nodelete");
2387  nodes.Clear("nodelete");
2388 }
2389 
2390 
2391 ////////////////////////////////////////////////////////////////////////////////
2392 /// Wrapper of all exporting methods
2393 /// Creates blank GDML file and fills it with gGeoManager structure converted
2394 /// to GDML structure of xml nodes
2395 
2396 void TGDMLWrite::WriteGDMLfile(TGeoManager * geomanager,
2397  TGeoVolume* volume,
2398  TList* materialsLst,
2399  const char* filename,
2400  TString option)
2401 {
2402  //option processing
2403  option.ToLower();
2404  if (option.Contains("g")) {
2405  SetG4Compatibility(kTRUE);
2406  Info("WriteGDMLfile", "Geant4 compatibility mode set");
2407  } else {
2408  SetG4Compatibility(kFALSE);
2409  }
2410  if (option.Contains("f")) {
2411  SetNamingSpeed(kfastButUglySufix);
2412  Info("WriteGDMLfile", "Fast naming convention with pointer suffix set");
2413  } else if (option.Contains("n")) {
2414  SetNamingSpeed(kwithoutSufixNotUniq);
2415  Info("WriteGDMLfile", "Naming without prefix set - be careful uniqness of name is not ensured");
2416  } else {
2417  SetNamingSpeed(kelegantButSlow);
2418  Info("WriteGDMLfile", "Potentially slow with incremental suffix naming convention set");
2419  }
2420 
2421  //local variables
2422  Int_t outputLayout = 1;
2423  const char * krootNodeName = "gdml";
2424  const char * knsRefGeneral = "http://www.w3.org/2001/XMLSchema-instance";
2425  const char * knsNameGeneral = "xsi";
2426  const char * knsRefGdml = "http://service-spi.web.cern.ch/service-spi/app/releases/GDML/schema/gdml.xsd";
2427  const char * knsNameGdml = "xsi:noNamespaceSchemaLocation";
2428 
2429  // First create engine
2430  fGdmlE = new TXMLEngine;
2431  fGdmlE->SetSkipComments(kTRUE);
2432 
2433  //create blank GDML file
2434  fGdmlFile = fGdmlE->NewDoc();
2435 
2436  //create root node and add it to blank GDML file
2437  XMLNodePointer_t rootNode = fGdmlE->NewChild(nullptr, nullptr, krootNodeName, nullptr);
2438  fGdmlE->DocSetRootElement(fGdmlFile, rootNode);
2439 
2440  //add namespaces to root node
2441  fGdmlE->NewNS(rootNode, knsRefGeneral, knsNameGeneral);
2442  fGdmlE->NewAttr(rootNode, nullptr, knsNameGdml, knsRefGdml);
2443 
2444  //initialize general lists and <define>, <solids>, <structure> nodes
2445  fIsotopeList = new StructLst;
2446  fElementList = new StructLst;
2447 
2448  fNameList = new NameLst;
2449 
2450  fDefineNode = fGdmlE->NewChild(nullptr, nullptr, "define", nullptr);
2451  fSolidsNode = fGdmlE->NewChild(nullptr, nullptr, "solids", nullptr);
2452  fStructureNode = fGdmlE->NewChild(nullptr, nullptr, "structure", nullptr);
2453  //========================
2454 
2455  //initialize list of accepted patterns for divisions (in ExtractVolumes)
2456  fAccPatt = new StructLst;
2457  fAccPatt->fLst["TGeoPatternX"] = kTRUE;
2458  fAccPatt->fLst["TGeoPatternY"] = kTRUE;
2459  fAccPatt->fLst["TGeoPatternZ"] = kTRUE;
2460  fAccPatt->fLst["TGeoPatternCylR"] = kTRUE;
2461  fAccPatt->fLst["TGeoPatternCylPhi"] = kTRUE;
2462  //========================
2463 
2464  //initialize list of rejected shapes for divisions (in ExtractVolumes)
2465  fRejShape = new StructLst;
2466  //this shapes are rejected because, it is not possible to divide trd2
2467  //in Y axis and while only trd2 object is imported from GDML
2468  //it causes a problem when TGeoTrd1 is divided in Y axis
2469  fRejShape->fLst["TGeoTrd1"] = kTRUE;
2470  fRejShape->fLst["TGeoTrd2"] = kTRUE;
2471  //=========================
2472 
2473  //Initialize global counters
2474  fActNameErr = 0;
2475  fVolCnt = 0;
2476  fPhysVolCnt = 0;
2477  fSolCnt = 0;
2478 
2479  //calling main extraction functions (with measuring time)
2480  time_t startT, endT;
2481  startT = time(nullptr);
2482  ExtractMatrices(geomanager->GetListOfGDMLMatrices());
2483  ExtractConstants(geomanager);
2484  fMaterialsNode = ExtractMaterials(materialsLst);
2485 
2486  Info("WriteGDMLfile", "Extracting volumes");
2487  ExtractVolumes(volume);
2488  Info("WriteGDMLfile", "%i solids added", fSolCnt);
2489  Info("WriteGDMLfile", "%i volumes added", fVolCnt);
2490  Info("WriteGDMLfile", "%i physvolumes added", fPhysVolCnt);
2491  ExtractSkinSurfaces(geomanager->GetListOfSkinSurfaces());
2492  ExtractBorderSurfaces(geomanager->GetListOfBorderSurfaces());
2493  ExtractOpticalSurfaces(geomanager->GetListOfOpticalSurfaces());
2494  endT = time(nullptr);
2495  //<gdml>
2496  fGdmlE->AddChild(rootNode, fDefineNode); // <define>...</define>
2497  fGdmlE->AddChild(rootNode, fMaterialsNode); // <materials>...</materials>
2498  fGdmlE->AddChild(rootNode, fSolidsNode); // <solids>...</solids>
2499  fGdmlE->AddChild(rootNode, fStructureNode); // <structure>...</structure>
2500  fGdmlE->AddChild(rootNode, CreateSetupN(fTopVolumeName.Data())); // <setup>...</setup>
2501  //</gdml>
2502  Double_t tdiffI = difftime(endT, startT);
2503  TString tdiffS = (tdiffI == 0 ? TString("< 1 s") : TString::Format("%.0lf s", tdiffI));
2504  Info("WriteGDMLfile", "Exporting time: %s", tdiffS.Data());
2505  //=========================
2506 
2507  //Saving document
2508  fGdmlE->SaveDoc(fGdmlFile, filename, outputLayout);
2509  Info("WriteGDMLfile", "File %s saved", filename);
2510  //cleaning
2511  fGdmlE->FreeDoc(fGdmlFile);
2512  //unset processing bits:
2513  UnsetTemporaryBits(geomanager);
2514  delete fGdmlE;
2515 }
2516 
2517 
2518 
2519 
2520 ////////////////////////////////////////////////////////////////////////////////
2521 /// Method extracting geometry structure recursively
2522 
2523 void TGDMLWrite::ExtractVolumes(TGeoVolume* volume)
2524 {
2525  XMLNodePointer_t volumeN, childN;
2526  TString volname, matname, solname, pattClsName, nodeVolNameBak;
2527  TGeoPatternFinder *pattFinder = nullptr;
2528  Bool_t isPattern = kFALSE;
2529  const TString fltPrecision = TString::Format("%%.%dg", fFltPrecision);
2530 
2531  //create the name for volume/assembly
2532  if (volume->IsTopVolume()) {
2533  //not needed a special function for generating name
2534  volname = volume->GetName();
2535  fTopVolumeName = volname;
2536  //register name to the pointer
2537  fNameList->fLst[TString::Format("%p", volume)] = volname;
2538  } else {
2539  volname = GenName(volume->GetName(), TString::Format("%p", volume));
2540  }
2541 
2542  //start to create main volume/assembly node
2543  if (volume->IsAssembly()) {
2544  volumeN = StartAssemblyN(volname);
2545  } else {
2546  //get reference material and add solid to <solids> + get name
2547  matname = fNameList->fLst[TString::Format("%p", volume->GetMaterial())];
2548  solname = ExtractSolid(volume->GetShape());
2549  //If solid is not supported or corrupted
2550  if (solname == "-1") {
2551  Info("ExtractVolumes", "ERROR! %s volume was not added, because solid is either not supported or corrupted",
2552  volname.Data());
2553  //set volume as missing volume
2554  fNameList->fLst[TString::Format("%p", volume)] = "missing_" + volname;
2555  return;
2556  }
2557  volumeN = StartVolumeN(volname, solname, matname);
2558 
2559  //divisionvol can't be in assembly
2560  pattFinder = volume->GetFinder();
2561  //if found pattern
2562  if (pattFinder) {
2563  pattClsName = TString::Format("%s", pattFinder->ClassName());
2564  TString shapeCls = TString::Format("%s", volume->GetShape()->ClassName());
2565  //if pattern in accepted pattern list and not in shape rejected list
2566  if ((fAccPatt->fLst[pattClsName] == kTRUE) &&
2567  (fRejShape->fLst[shapeCls] != kTRUE)) {
2568  isPattern = kTRUE;
2569  }
2570  }
2571  }
2572  //get all nodes in volume
2573  TObjArray *nodeLst = volume->GetNodes();
2574  TIter next(nodeLst);
2575  TGeoNode *geoNode;
2576  Int_t nCnt = 0;
2577  //loop through all nodes
2578  while ((geoNode = (TGeoNode *) next())) {
2579  //get volume of current node and if not processed then process it
2580  TGeoVolume * subvol = geoNode->GetVolume();
2581  if (subvol->TestAttBit(fgkProcBitVol) == kFALSE) {
2582  subvol->SetAttBit(fgkProcBitVol);
2583  ExtractVolumes(subvol);
2584  }
2585 
2586  //volume of this node has to exist because it was processed recursively
2587  TString nodevolname = fNameList->fLst[TString::Format("%p", geoNode->GetVolume())];
2588  if (nodevolname.Contains("missing_")) {
2589  continue;
2590  }
2591  if (nCnt == 0) { //save name of the first node for divisionvol
2592  nodeVolNameBak = nodevolname;
2593  }
2594 
2595  if (isPattern == kFALSE) {
2596  //create name for node
2597  TString nodename, posname, rotname;
2598  nodename = GenName(geoNode->GetName(), TString::Format("%p", geoNode));
2599  nodename = nodename + "in" + volname;
2600 
2601  //create name for position and clear rotation
2602  posname = nodename + "pos";
2603  rotname = "";
2604 
2605  //position
2606  const Double_t * pos = geoNode->GetMatrix()->GetTranslation();
2607  Xyz nodPos;
2608  nodPos.x = pos[0];
2609  nodPos.y = pos[1];
2610  nodPos.z = pos[2];
2611  childN = CreatePositionN(posname.Data(), nodPos);
2612  fGdmlE->AddChild(fDefineNode, childN); //adding node to <define> node
2613  //Deal with reflection
2614  XMLNodePointer_t scaleN = nullptr;
2615  Double_t lx, ly, lz;
2616  Double_t xangle = 0;
2617  Double_t zangle = 0;
2618  lx = geoNode->GetMatrix()->GetRotationMatrix()[0];
2619  ly = geoNode->GetMatrix()->GetRotationMatrix()[4];
2620  lz = geoNode->GetMatrix()->GetRotationMatrix()[8];
2621  if (geoNode->GetMatrix()->IsReflection()
2622  && TMath::Abs(lx) == 1 && TMath::Abs(ly) == 1 && TMath::Abs(lz) == 1) {
2623  scaleN = fGdmlE->NewChild(nullptr, nullptr, "scale", nullptr);
2624  fGdmlE->NewAttr(scaleN, nullptr, "name", (nodename + "scl").Data());
2625  fGdmlE->NewAttr(scaleN, nullptr, "x", TString::Format(fltPrecision.Data(), lx));
2626  fGdmlE->NewAttr(scaleN, nullptr, "y", TString::Format(fltPrecision.Data(), ly));
2627  fGdmlE->NewAttr(scaleN, nullptr, "z", TString::Format(fltPrecision.Data(), lz));
2628  //experimentally found out, that rotation should be updated like this
2629  if (lx == -1) {
2630  zangle = 180;
2631  }
2632  if (lz == -1) {
2633  xangle = 180;
2634  }
2635  }
2636 
2637  //rotation
2638  TGDMLWrite::Xyz lxyz = GetXYZangles(geoNode->GetMatrix()->GetRotationMatrix());
2639  lxyz.x -= xangle;
2640  lxyz.z -= zangle;
2641  if ((lxyz.x != 0.0) || (lxyz.y != 0.0) || (lxyz.z != 0.0)) {
2642  rotname = nodename + "rot";
2643  childN = CreateRotationN(rotname.Data(), lxyz);
2644  fGdmlE->AddChild(fDefineNode, childN); //adding node to <define> node
2645  }
2646 
2647  //create physvol for main volume/assembly node
2648  childN = CreatePhysVolN(geoNode->GetName(), geoNode->GetNumber(), nodevolname.Data(), posname.Data(), rotname.Data(), scaleN);
2649  fGdmlE->AddChild(volumeN, childN);
2650  }
2651  nCnt++;
2652  }
2653  //create only one divisionvol node
2654  if (isPattern && pattFinder) {
2655  //retrieve attributes of division
2656  Int_t ndiv, divaxis;
2657  Double_t offset, width, xlo, xhi;
2658  TString axis, unit;
2659 
2660  ndiv = pattFinder->GetNdiv();
2661  width = pattFinder->GetStep();
2662 
2663  divaxis = pattFinder->GetDivAxis();
2664  volume->GetShape()->GetAxisRange(divaxis, xlo, xhi);
2665 
2666  //compute relative start (not positional)
2667  offset = pattFinder->GetStart() - xlo;
2668  axis = GetPattAxis(divaxis, pattClsName, unit);
2669 
2670  //create division node
2671  childN = CreateDivisionN(offset, width, ndiv, axis.Data(), unit.Data(), nodeVolNameBak.Data());
2672  fGdmlE->AddChild(volumeN, childN);
2673  }
2674 
2675  fVolCnt++;
2676  //add volume/assembly node into the <structure> node
2677  fGdmlE->AddChild(fStructureNode, volumeN);
2678 }