Belle II Software  release-08-01-10
VRMLWriterModule.cc
1 /**************************************************************************
2  * basf2 (Belle II Analysis Software Framework) *
3  * Author: The Belle II Collaboration *
4  * *
5  * See git log for contributors and copyright holders. *
6  * This file is licensed under LGPL-3.0, see LICENSE.md. *
7  **************************************************************************/
8 
9 #include <geometry/modules/vrmlWriter/VRMLWriterModule.h>
10 #include <geometry/GeometryManager.h>
11 
12 #include "G4PhysicalVolumeStore.hh"
13 #include "G4LogicalVolumeStore.hh"
14 #include "G4SolidStore.hh"
15 #include "G4VPhysicalVolume.hh"
16 #include "G4LogicalVolume.hh"
17 #include "G4VSolid.hh"
18 #include "G4VisAttributes.hh"
19 #include "G4VisExtent.hh"
20 #include "G4Material.hh"
21 #include "G4ThreeVector.hh"
22 #include "G4RotationMatrix.hh"
23 #include "G4Transform3D.hh"
24 #include "G4VPVParameterisation.hh"
25 #include <G4Tubs.hh>
26 #include <G4Box.hh>
27 #include <G4Polyhedron.hh>
28 #include <G4DisplacedSolid.hh>
29 
30 #include <iomanip>
31 
32 using namespace Belle2;
33 
34 REG_MODULE(VRMLWriter);
35 
36 VRMLWriterModule::VRMLWriterModule()
37 {
38  //Set module properties and the description
39  setDescription("Write the detector geometry in a hierarchical VRML format.");
40 
41  //Parameter definition
42  addParam("outputFile", m_Filename, "Output filename", std::string("belle2.wrl"));
43 }
44 
45 void VRMLWriterModule::initialize()
46 {
47  m_First = true;
48 }
49 
50 void VRMLWriterModule::event()
51 {
52  if (!m_First) return;
53  m_First = false;
54  G4VPhysicalVolume* topVol = geometry::GeometryManager::getInstance().getTopVolume();
55  if (!topVol) {
56  B2ERROR("No geometry found: add the Geometry module to the path before the VRMLWriter module.");
57  return;
58  }
59  m_File.open(m_Filename, std::ios_base::trunc);
60  writePreamble();
61 
62  G4PhysicalVolumeStore* pvStore = G4PhysicalVolumeStore::GetInstance();
63  G4LogicalVolumeStore* lvStore = G4LogicalVolumeStore::GetInstance();
64  G4SolidStore* solidStore = G4SolidStore::GetInstance();
65 
66  m_PVName = new std::vector<std::string>(pvStore->size(), "");
67  m_LVName = new std::vector<std::string>(lvStore->size(), "");
68  m_SolidName = new std::vector<std::string>(solidStore->size(), "");
69 
70  // Assign legal and unique names to each used physical volume, logical volume and solid
71  for (G4VPhysicalVolume* physVol : *pvStore) {
72  int pvIndex = std::find(pvStore->begin(), pvStore->end(), physVol) - pvStore->begin();
73  if ((*m_PVName)[pvIndex].length() == 0)
74  assignName(m_PVName, pvIndex, (*pvStore)[pvIndex]->GetName(), 0);
75  G4LogicalVolume* logVol = physVol->GetLogicalVolume();
76  int lvIndex = std::find(lvStore->begin(), lvStore->end(), logVol) - lvStore->begin();
77  if ((*m_LVName)[lvIndex].length() == 0)
78  assignName(m_LVName, lvIndex, (*lvStore)[lvIndex]->GetName(), 1);
79  G4VSolid* solid = logVol->GetSolid();
80  int solidIndex = std::find(solidStore->begin(), solidStore->end(), solid) - solidStore->begin();
81  if ((*m_SolidName)[solidIndex].length() == 0)
82  assignName(m_SolidName, solidIndex, (*solidStore)[solidIndex]->GetName(), 2);
83  }
84 
85  // Write all explicitly-referenced solids as PROTOs (replicas are written later)
86  // Use implicit prefix "solid_" to avoid name clashes with logical- and physical-volume names
87  m_IsCylinder = new std::vector<bool>(solidStore->size(), false);
88  for (unsigned int solidIndex = 0; solidIndex < solidStore->size(); ++solidIndex) {
89  if ((*m_SolidName)[solidIndex].length() != 0) {
90  if ((*solidStore)[solidIndex]->GetEntityType() == "G4Tubs") {
91  auto* tube = (G4Tubs*)((*solidStore)[solidIndex]);
92  (*m_IsCylinder)[solidIndex] = ((tube->GetInnerRadius() == 0.0) && (tube->GetDeltaPhiAngle() == 2.0 * M_PI));
93  }
94  describeSolid((*solidStore)[solidIndex], (*m_SolidName)[solidIndex], (*m_IsCylinder)[solidIndex]);
95  }
96  }
97 
98  // Recursively write all physical volumes (as DEFs) and logical volumes (as PROTOs).
99  // Deepest volumes are written first; top volume is written last. Use implicit prefix
100  // "lv_" for logical-volume names to avoid name clashes with solid and physical-volume names.
101  m_PVIndex = new std::vector<std::vector<int> >(lvStore->size(), std::vector<int>());
102  m_LVWritten = new std::vector<bool>(lvStore->size(), false);
103  m_PVWritten = new std::vector<bool>(pvStore->size(), false);
104  describePhysicalVolume(topVol);
105 
106  // Now tell VRML to draw the top physical volume (and, recursively, all daughters)
107  int pvIndex = std::find(pvStore->begin(), pvStore->end(), topVol) - pvStore->begin();
108  int lvIndex = std::find(lvStore->begin(), lvStore->end(), topVol->GetLogicalVolume()) - lvStore->begin();
109  m_File << "# PhysicalVolume " << topVol->GetName() << std::endl <<
110  "DEF " << (*m_PVName)[pvIndex] << " Transform {" << std::endl <<
111  " children lv_" << (*m_LVName)[lvIndex] << "{}" << std::endl <<
112  "}" << std::endl;
113  m_File.close();
114  B2INFO("VRML written to " << m_Filename);
115 
116  delete m_PVName;
117  delete m_LVName;
118  delete m_SolidName;
119  delete m_IsCylinder;
120  delete m_PVIndex;
121  delete m_LVWritten;
122  delete m_PVWritten;
123 
124 }
125 
126 void VRMLWriterModule::assignName(std::vector<std::string>* names, unsigned int index, const G4String& originalName, int select)
127 {
128  G4SolidStore* solidStore = G4SolidStore::GetInstance();
129  G4LogicalVolumeStore* lvStore = G4LogicalVolumeStore::GetInstance();
130  G4PhysicalVolumeStore* pvStore = G4PhysicalVolumeStore::GetInstance();
131 
132  G4String name = originalName;
133  if (name.length() == 0) { name = "anonymous"; }
134  // Replace problematic characters with underscore (there may be more!)
135  for (char c : " .,:;?'\"*+-=|^!/@#$\\%{}[]()<>") std::replace(name.begin(), name.end(), c, '_');
136  // Avoid duplicate names for solids that will be written to VRML file
137  for (int j = (int)index - 1; j >= 0; --j) {
138  if ((*names)[j].length() == 0) continue;
139  int match = 0;
140  switch (select) {
141  case 0:
142  match = (*pvStore)[j]->GetName().compare((*pvStore)[index]->GetName()); break;
143  case 1:
144  match = (*lvStore)[j]->GetName().compare((*lvStore)[index]->GetName()); break;
145  case 2:
146  match = (*solidStore)[j]->GetName().compare((*solidStore)[index]->GetName()); break;
147  }
148  if (match == 0) {
149  if (name.length() == (*names)[j].length()) {
150  (*names)[j].append("_1");
151  }
152  int n = std::stoi((*names)[j].substr(name.length() + 1), nullptr);
153  name.append("_");
154  name.append(std::to_string(n + 1));
155  break;
156  }
157  }
158  (*names)[index] = name;
159 }
160 
161 void VRMLWriterModule::describeSolid(G4VSolid* solid, const std::string& solidName, bool isCylinder)
162 {
163  m_File << "# Solid " << solid->GetName() << " of type " << solid->GetEntityType() << std::endl;
164  if (isCylinder) {
165  auto* tube = (G4Tubs*)solid;
166  // VRML cylinder is along y axis but G4Tubs is along z axis => rotate in logical volume
167  m_File << "PROTO solid_" << solidName << " [ ] {" << std::endl <<
168  " Cylinder {" << std::endl << std::setprecision(10) <<
169  " radius " << tube->GetOuterRadius() << std::endl <<
170  " height " << tube->GetZHalfLength() * 2.0 << std::endl <<
171  " }" << std::endl <<
172  "}" << std::endl;
173  } else if (solid->GetEntityType() == "G4Box") {
174  auto* box = (G4Box*)solid;
175  m_File << "PROTO solid_" << solidName << " [ ] {" << std::endl <<
176  " Box {" << std::endl << std::setprecision(10) <<
177  " size " << box->GetXHalfLength() * 2.0 << " " <<
178  box->GetYHalfLength() * 2.0 << " " <<
179  box->GetZHalfLength() * 2.0 << std::endl <<
180  " }" << std::endl <<
181  "}" << std::endl;
182  } else if ((solid->GetEntityType() == "G4IntersectionSolid") ||
183  (solid->GetEntityType() == "G4UnionSolid") ||
184  (solid->GetEntityType() == "G4SubtractionSolid") ||
185  (solid->GetEntityType() == "G4BooleanSolid")) {
186  HepPolyhedron* polyhedron = getBooleanSolidPolyhedron(solid);
187  auto* g4polyhedron = new G4Polyhedron(*polyhedron);
188  writePolyhedron(g4polyhedron, solidName);
189  delete polyhedron;
190  delete g4polyhedron;
191  } else {
192  writePolyhedron(solid->GetPolyhedron(), solidName);
193  }
194 }
195 
196 void VRMLWriterModule::describeLogicalVolume(G4LogicalVolume* logVol, const std::string& lvName, const std::string& solidName,
197  bool isCylinder)
198 {
199  G4LogicalVolumeStore* lvStore = G4LogicalVolumeStore::GetInstance();
200  int lvIndex = std::find(lvStore->begin(), lvStore->end(), logVol) - lvStore->begin();
201  if ((*m_LVWritten)[lvIndex]) return; // at most one PROTO in VRML for each logical volume
202  (*m_LVWritten)[lvIndex] = true;
203  G4Color color(0.0, 1.0, 0.0, 0.5); // default is semi-transparent green
204  if ((lvName.compare(0, 23, "eclBarrelCrystalLogical") == 0) ||
205  (lvName.compare(0, 20, "eclFwdCrystalLogical") == 0) ||
206  (lvName.compare(0, 20, "eclBwdCrystalLogical") == 0)) {
207  color = G4Color(1.0, 0.25, 0.0, 0.7); // orange since ECL crystals have no G4VisAttribute :(
208  }
209  std::string visible = "";
210  G4String materialName = logVol->GetMaterial()->GetName();
211  // Hide containers that have vacuum, air or gas
212  if (materialName == "Vacuum") visible = "#";
213  if (materialName == "G4_AIR") visible = "#";
214  if (materialName == "CDCGas") visible = "#";
215  if (materialName == "ColdAir") visible = "#";
216  if (materialName == "STR-DryAir") visible = "#";
217  if (materialName == "TOPAir") visible = "#";
218  if (materialName == "TOPVacuum") visible = "#";
219  const G4VisAttributes* visAttr = logVol->GetVisAttributes();
220  if (visAttr) {
221  color = const_cast<G4Color&>(logVol->GetVisAttributes()->GetColor());
222  if (!(visAttr->IsVisible())) visible = "#";
223  } else {
224  visible = "#";
225  }
226  bool hasDaughters = (visible.length() == 0) || ((*m_PVIndex)[lvIndex].size() > 0);
227  if (logVol->GetSensitiveDetector() != nullptr) visible = "";
228  m_File << "# LogicalVolume " << logVol->GetName() << " containing " << materialName << std::endl <<
229  "PROTO lv_" << lvName << " [ ] {" << std::endl <<
230  " Group {" << std::endl <<
231  " " << (hasDaughters ? "" : "#") << "children [" << std::endl;
232  if (isCylinder) {
233  // VRML cylinder is along y axis but G4Tubs is along z axis => rotate here
234  m_File << " " << visible << "Transform {" << std::endl <<
235  " " << visible << " rotation 1 0 0 " << M_PI_2 << std::endl <<
236  " " << visible << " children [" << std::endl;
237  visible.append(" ");
238  }
239  m_File << " " << visible << "Shape {" << std::endl <<
240  " " << visible << " appearance Appearance {" << std::endl <<
241  " " << visible << " material Material {" << std::endl <<
242  " " << visible << " diffuseColor " << color.GetRed() << " " <<
243  color.GetGreen() << " " <<
244  color.GetBlue() << std::endl <<
245  " " << visible << " " << (color.GetAlpha() == 1.0 ? "#" : "") <<
246  "transparency " << 1.0 - color.GetAlpha() << std::endl <<
247  " " << visible << " }" << std::endl <<
248  " " << visible << " }" << std::endl <<
249  " " << visible << " geometry solid_" << solidName << "{}" << std::endl <<
250  " " << visible << "}" << std::endl;
251  if (isCylinder) {
252  visible.resize(visible.length() - 2);
253  m_File << " " << visible << " ]" << std::endl <<
254  " " << visible << "}" << std::endl;
255  }
256  G4PhysicalVolumeStore* pvStore = G4PhysicalVolumeStore::GetInstance();
257  for (int daughter : (*m_PVIndex)[lvIndex]) {
258  G4VPhysicalVolume* physVol = (*pvStore)[daughter];
259  std::string pvNameDaughter = (*m_PVName)[daughter];
260  int lvDaughter = std::find(lvStore->begin(), lvStore->end(), physVol->GetLogicalVolume()) - lvStore->begin();
261  std::string lvNameDaughter = (*m_LVName)[lvDaughter];
262  if (physVol->IsReplicated()) {
263  G4VSolid* solid = logVol->GetSolid();
264  EAxis axis;
265  G4int nReplicas;
266  G4double width;
267  G4double offset;
268  G4bool consuming;
269  physVol->GetReplicationData(axis, nReplicas, width, offset, consuming);
270  G4VPVParameterisation* physParameterisation = physVol->GetParameterisation();
271  if (physParameterisation) { // parameterised volume
272  for (int replica = 0; replica < nReplicas; ++replica) {
273  std::string pvNameReplica = pvNameDaughter;
274  std::string lvNameReplica = lvNameDaughter;
275  pvNameReplica.append("_");
276  pvNameReplica.append(std::to_string(replica));
277  physVol->SetCopyNo(replica);
278  physParameterisation->ComputeTransformation(replica, physVol);
279  G4VSolid* solidReplica = physParameterisation->ComputeSolid(replica, physVol);
280  if (solidReplica != solid) { // not sure if this works ...
281  lvNameReplica.append("_");
282  lvNameReplica.append(std::to_string(replica));
283  }
284  writePhysicalVolume(physVol, pvNameReplica, lvNameReplica, false);
285  }
286  } else { // plain replicated volume
287  G4ThreeVector originalTranslation = physVol->GetTranslation();
288  G4RotationMatrix* originalRotation = physVol->GetRotation();
289  for (int replica = 0; replica < nReplicas; ++replica) {
290  G4ThreeVector translation; // No translation
291  G4RotationMatrix rotation; // No rotation
292  std::string pvNameReplica = pvNameDaughter;
293  std::string lvNameReplica = lvNameDaughter;
294  pvNameReplica.append("_");
295  pvNameReplica.append(std::to_string(replica));
296  physVol->SetCopyNo(replica);
297  switch (axis) {
298  default:
299  case kXAxis:
300  translation.setX(width * (replica - 0.5 * (nReplicas - 1)));
301  physVol->SetTranslation(translation);
302  break;
303  case kYAxis:
304  translation.setY(width * (replica - 0.5 * (nReplicas - 1)));
305  physVol->SetTranslation(translation);
306  break;
307  case kZAxis:
308  translation.setZ(width * (replica - 0.5 * (nReplicas - 1)));
309  physVol->SetTranslation(translation);
310  break;
311  case kRho:
312  if (solid->GetEntityType() == "G4Tubs") {
313  lvNameReplica.append("_");
314  lvNameReplica.append(std::to_string(replica));
315  } else {
316  B2WARNING("Built-in volumes replicated along radius for " << solid->GetEntityType() <<
317  " (solid " << solid->GetName() << ") are not visualisable.");
318  }
319  break;
320  case kPhi:
321  physVol->SetRotation(&(rotation.rotateZ(-(offset + (replica + 0.5) * width))));
322  break;
323  }
324  writePhysicalVolume(physVol, pvNameReplica, lvNameReplica, false);
325  }
326  // Restore originals...
327  physVol->SetTranslation(originalTranslation);
328  physVol->SetRotation(originalRotation);
329  }
330  } else {
331  writePhysicalVolume(physVol, pvNameDaughter, lvNameDaughter, (*m_PVWritten)[daughter]);
332  }
333  (*m_PVWritten)[daughter] = true;
334  }
335  G4VisExtent extent = logVol->GetSolid()->GetExtent();
336  double xMin = extent.GetXmin();
337  double xMax = extent.GetXmax();
338  double yMin = extent.GetYmin();
339  double yMax = extent.GetYmax();
340  double zMin = extent.GetZmin();
341  double zMax = extent.GetZmax();
342  double xCenter = 0.5 * (xMin + xMax);
343  double yCenter = 0.5 * (yMin + yMax);
344  double zCenter = 0.5 * (zMin + zMax);
345  bool atOrigin = (std::fabs(xCenter) + std::fabs(yCenter) + std::fabs(zCenter) > 1.0e-12);
346  m_File << " " << (hasDaughters ? "" : "#") << "]" << std::endl <<
347  " " << (atOrigin ? "" : "#") <<
348  "bboxCenter " << xCenter << " " << yCenter << " " << zCenter << std::endl <<
349  " bboxSize " << xMax - xMin << " " <<
350  (isCylinder ? zMax - zMin : yMax - yMin) << " " <<
351  (isCylinder ? yMax - yMin : zMax - zMin) << std::endl <<
352  " }" << std::endl <<
353  "} # end of " << lvName << std::endl;
354 }
355 
356 void VRMLWriterModule::describePhysicalVolume(G4VPhysicalVolume* physVol)
357 {
358  G4LogicalVolume* logVol = physVol->GetLogicalVolume();
359  G4VSolid* solid = logVol->GetSolid();
360  G4LogicalVolumeStore* lvStore = G4LogicalVolumeStore::GetInstance();
361  int lvIndex = std::find(lvStore->begin(), lvStore->end(), logVol) - lvStore->begin();
362  std::string lvName = (*m_LVName)[lvIndex];
363  if (physVol->IsReplicated()) {
364  G4SolidStore* solidStore = G4SolidStore::GetInstance();
365  int solidIndex = std::find(solidStore->begin(), solidStore->end(), solid) - solidStore->begin();
366  std::string solidName = (*m_SolidName)[solidIndex];
367  EAxis axis;
368  G4int nReplicas;
369  G4double width;
370  G4double offset;
371  G4bool consuming;
372  physVol->GetReplicationData(axis, nReplicas, width, offset, consuming);
373  G4VPVParameterisation* physParameterisation = physVol->GetParameterisation();
374  if (physParameterisation) { // user-parameterised volume
375  for (int replica = 0; replica < nReplicas; ++replica) {
376  std::string lvNameReplica = lvName;
377  G4VSolid* solidReplica = physParameterisation->ComputeSolid(replica, physVol);
378  if (solidReplica != solid) { // not sure if this works ...
379  solidReplica->ComputeDimensions(physParameterisation, replica, physVol);
380  std::string solidNameReplica = solidName;
381  solidNameReplica.append("_");
382  solidNameReplica.append(std::to_string(replica));
383  lvNameReplica.append("_");
384  lvNameReplica.append(std::to_string(replica));
385  describeSolid(solidReplica, solidNameReplica, (*m_IsCylinder)[solidIndex]);
386  describeLogicalVolume(logVol, lvNameReplica, solidNameReplica, (*m_IsCylinder)[solidIndex]);
387  }
388  descendAndDescribe(physVol, lvNameReplica, replica);
389  }
390  } else { // plain replicated volume - "divisions"
391  for (int replica = 0; replica < nReplicas; ++replica) {
392  if (axis == kRho) {
393  if (solid->GetEntityType() == "G4Tubs") {
394  double originalRMin = ((G4Tubs*)solid)->GetInnerRadius();
395  double originalRMax = ((G4Tubs*)solid)->GetOuterRadius();
396  ((G4Tubs*)solid)->SetInnerRadius(offset + width * replica);
397  ((G4Tubs*)solid)->SetOuterRadius(offset + width * (replica + 1));
398  std::string solidNameReplica = solidName;
399  solidNameReplica.append("_");
400  solidNameReplica.append(std::to_string(replica));
401  std::string lvNameReplica = lvName;
402  lvNameReplica.append("_");
403  lvNameReplica.append(std::to_string(replica));
404  describeSolid(solid, solidNameReplica, (*m_IsCylinder)[solidIndex]);
405  ((G4Tubs*)solid)->SetInnerRadius(originalRMin);
406  ((G4Tubs*)solid)->SetOuterRadius(originalRMax);
407  describeLogicalVolume(logVol, lvNameReplica, solidNameReplica, (*m_IsCylinder)[solidIndex]);
408  descendAndDescribe(physVol, lvNameReplica, replica);
409  } else {
410  B2WARNING("Built-in volumes replicated along radius for " << solid->GetEntityType() <<
411  " (solid " << solid->GetName() << ") are not visualisable.");
412  }
413  } else {
414  descendAndDescribe(physVol, lvName, replica);
415  }
416  }
417  }
418  } else { // non-replicated volume
419  descendAndDescribe(physVol, lvName, -1);
420  }
421 }
422 
423 void VRMLWriterModule::descendAndDescribe(G4VPhysicalVolume* physVol, const std::string& lvName, int replica)
424 {
425  G4PhysicalVolumeStore* pvStore = G4PhysicalVolumeStore::GetInstance();
426  G4LogicalVolumeStore* lvStore = G4LogicalVolumeStore::GetInstance();
427  G4SolidStore* solidStore = G4SolidStore::GetInstance();
428 
429  int pvIndex = std::find(pvStore->begin(), pvStore->end(), physVol) - pvStore->begin();
430  if ((*m_PVWritten)[pvIndex]) return;
431 
432  // Descend to the leaves of the tree
433  G4LogicalVolume* logVol = physVol->GetLogicalVolume();
434  int lvIndex = std::find(lvStore->begin(), lvStore->end(), logVol) - lvStore->begin();
435  for (size_t daughter = 0; daughter < logVol->GetNoDaughters(); ++daughter) {
436  G4VPhysicalVolume* dPhysVol = logVol->GetDaughter(daughter);
437  (*m_PVIndex)[lvIndex].push_back(std::find(pvStore->begin(), pvStore->end(), dPhysVol) - pvStore->begin());
438  describePhysicalVolume(dPhysVol);
439  }
440  // Write out the physical volume and its corresponding logical volume as we ascend the recursive tree
441  int solidIndex = std::find(solidStore->begin(), solidStore->end(), logVol->GetSolid()) - solidStore->begin();
442  if (replica <= 0) describeLogicalVolume(logVol, lvName, (*m_SolidName)[solidIndex], (*m_IsCylinder)[solidIndex]);
443 }
444 
445 void VRMLWriterModule::writePreamble()
446 {
447  m_File << "#VRML V2.0 utf8" << std::endl << std::endl <<
448  "WorldInfo {" << std::endl <<
449  " info [ \"(c) The Belle II Collaboration\" ]" << std::endl <<
450  " title \"The Belle II Detector\"" << std::endl <<
451  "}" << std::endl <<
452  "Viewpoint {" << std::endl <<
453  " position 0 0 15000" << std::endl <<
454  " description \"front\"" << std::endl <<
455  "}" << std::endl <<
456  "Viewpoint {" << std::endl <<
457  " position 0 0 -15000" << std::endl <<
458  " orientation 0 1 0 3.141593" << std::endl <<
459  " description \"back\"" << std::endl <<
460  "}" << std::endl <<
461  "Viewpoint {" << std::endl <<
462  " position 0 15000 0" << std::endl <<
463  " orientation 1 0 0 -1.570796" << std::endl <<
464  " description \"top\"" << std::endl <<
465  "}" << std::endl <<
466  "Viewpoint {" << std::endl <<
467  " position 0 -15000 0" << std::endl <<
468  " orientation 1 0 0 1.570796" << std::endl <<
469  " description \"bottom\"" << std::endl <<
470  "}" << std::endl <<
471  "Viewpoint {" << std::endl <<
472  " position 15000 0 0" << std::endl <<
473  " orientation 0 1 0 1.570796" << std::endl <<
474  " description \"right\"" << std::endl <<
475  "}" << std::endl <<
476  "Viewpoint {" << std::endl <<
477  " position -15000 0 0" << std::endl <<
478  " orientation 0 1 0 -1.570796" << std::endl <<
479  " description \"left\"" << std::endl <<
480  "}" << std::endl;
481 }
482 
483 void VRMLWriterModule::writePolyhedron(const G4Polyhedron* polyhedron, const std::string& name)
484 {
485  if (polyhedron) {
486  m_File << "PROTO solid_" << name << " [ ] {" << std::endl <<
487  " IndexedFaceSet {" << std::endl <<
488  " coord Coordinate {" << std::endl <<
489  " point [" << std::endl;
490  for (int j = 1; j <= polyhedron->GetNoVertices(); ++j) {
491  m_File << " " << polyhedron->GetVertex(j).x() << " " <<
492  polyhedron->GetVertex(j).y() << " " <<
493  polyhedron->GetVertex(j).z() << std::endl;
494  }
495  m_File << " ]" << std::endl <<
496  " }" << std::endl <<
497  " coordIndex [" << std::endl;
498  for (int j = 0; j < polyhedron->GetNoFacets(); ++j) {
499  G4bool notLastEdge = true;
500  G4int ndx = -1, edgeFlag = 1;
501  m_File << " ";
502  do {
503  notLastEdge = polyhedron->GetNextVertexIndex(ndx, edgeFlag);
504  m_File << " " << ndx - 1;
505  } while (notLastEdge);
506  m_File << " -1" << std::endl;
507  }
508  m_File << " ]" << std::endl <<
509  " }" << std::endl <<
510  "}" << std::endl;
511  } else {
512  B2INFO("Polyhedron representation of solid " << name << " cannot be created");
513  }
514 }
515 
516 void VRMLWriterModule::writePhysicalVolume(const G4VPhysicalVolume* physVol, const std::string& pvName, const std::string& lvName,
517  bool written)
518 {
519  if (written) {
520  m_File << " USE " << pvName << std::endl; // this never happens, as it turns out
521  } else {
522  m_File << " # PhysicalVolume " << physVol->GetName() << " copy " << physVol->GetCopyNo() << std::endl <<
523  " DEF " << pvName << " Transform {" << std::endl; // DEF provides object name when importing VRML into Unity
524  G4RotationMatrix* rot = physVol->GetObjectRotation();
525  G4ThreeVector move = physVol->GetObjectTranslation();
526  double yaw = std::atan2(rot->yx(), rot->xx());
527  double pitch = -std::asin(rot->zx());
528  double roll = std::atan2(rot->zy(), rot->zz());
529  std::cout << physVol->GetName() << " translation: " << move.x() << "," << move.y() << "," << move.z() << std::endl;
530  std::cout << physVol->GetName() << " rotation: " << yaw * 180.0 / M_PI << "," << pitch * 180.0 / M_PI << "," << roll * 180.0 / M_PI
531  << std::endl;
532  if (move.mag2() != 0.0) {
533  m_File << " translation " << std::setprecision(10) <<
534  move.x() << " " << move.y() << " " << move.z() << std::endl;
535  }
536  if (!(rot->isIdentity())) {
537  double trace = std::max(-1.0, std::min(3.0, rot->xx() + rot->yy() + rot->zz()));
538  double angle = std::acos(0.5 * (trace - 1.0));
539  G4ThreeVector rotA(rot->zy() - rot->yz(), rot->xz() - rot->zx(), rot->yx() - rot->xy());
540  if ((rotA.x() == 0.0) && (rotA.y() == 0.0) && (rotA.z() == 0.0)) {
541  // The assignment along x, y or z axis in Hep3Vector::axis() is wrong. The proper axis of
542  // rotation is given by kernel(r - I), which is the cross product of two non-proportional rows
543  if (rot->xx() > 0.0) {
544  rotA.setX((rot->yy() - 1.0) * (rot->zz() - 1.0) - rot->yz() * rot->zy());
545  rotA.setY(rot->yz() * rot->zx() - rot->yx() * (rot->zz() - 1.0));
546  rotA.setZ(rot->yx() * rot->zy() - rot->zx() * (rot->yy() - 1.0));
547  } else if (rot->yy() > 0.0) {
548  rotA.setX(rot->xy() * (rot->zz() - 1.0) - rot->xz() * rot->zy());
549  rotA.setY(rot->xz() * rot->zx() - (rot->xx() - 1.0) * (rot->zz() - 1.0));
550  rotA.setZ((rot->xx() - 1.0) * rot->zy() - rot->xy() * rot->zx());
551  } else {
552  rotA.setX(rot->xy() * rot->yz() - rot->xz() * (rot->yy() - 1.0));
553  rotA.setY(rot->xz() * rot->yx() - rot->yz() * (rot->xx() - 1.0));
554  rotA.setZ((rot->xx() - 1.0) * (rot->yy() - 1.0) - rot->xy() * rot->yx());
555  }
556  }
557  rotA.setMag(1.0);
558  m_File << " rotation " << std::setprecision(10) <<
559  rotA.x() << " " << rotA.y() << " " << rotA.z() << " " << angle << std::endl;
560  }
561  m_File << " children lv_" << lvName << "{}" << std::endl <<
562  " }" << std::endl;
563  }
564 }
565 
566 // The code in GEANT4 geometry/solids/Boolean/src/G4BooleanSolid.cc is buggy so avoid it.
567 // Recursively combine the polyhedrons of two solids to make a resulting polyhedron.
568 HepPolyhedron* VRMLWriterModule::getBooleanSolidPolyhedron(G4VSolid* solid)
569 {
570  G4VSolid* solidA = solid->GetConstituentSolid(0);
571  G4VSolid* solidB = solid->GetConstituentSolid(1);
572  HepPolyhedron* polyhedronA = nullptr;
573  if ((solidA->GetEntityType() == "G4IntersectionSolid") ||
574  (solidA->GetEntityType() == "G4UnionSolid") ||
575  (solidA->GetEntityType() == "G4SubtractionSolid") ||
576  (solidA->GetEntityType() == "G4BooleanSolid")) {
577  polyhedronA = getBooleanSolidPolyhedron(solidA);
578  } else {
579  polyhedronA = new HepPolyhedron(*(solidA->GetPolyhedron()));
580  }
581  HepPolyhedron* polyhedronB = nullptr;
582  G4VSolid* solidB2 = solidB;
583  if (solidB->GetEntityType() == "G4DisplacedSolid") {
584  solidB2 = ((G4DisplacedSolid*)solidB)->GetConstituentMovedSolid();
585  }
586  if ((solidB2->GetEntityType() == "G4IntersectionSolid") ||
587  (solidB2->GetEntityType() == "G4UnionSolid") ||
588  (solidB2->GetEntityType() == "G4SubtractionSolid") ||
589  (solidB2->GetEntityType() == "G4BooleanSolid")) {
590  polyhedronB = getBooleanSolidPolyhedron(solidB2);
591  if (solidB != solidB2) { // was solidB a G4DisplacedSolid?
592  polyhedronB->Transform(G4Transform3D(
593  ((G4DisplacedSolid*)solidB)->GetObjectRotation(),
594  ((G4DisplacedSolid*)solidB)->GetObjectTranslation()));
595  }
596  } else {
597  polyhedronB = new HepPolyhedron(*(solidB->GetPolyhedron()));
598  }
599  auto* result = new HepPolyhedron();
600  if (solid->GetEntityType() == "G4UnionSolid") {
601  *result = polyhedronA->add(*polyhedronB);
602  } else if (solid->GetEntityType() == "G4SubtractionSolid") {
603  *result = polyhedronA->subtract(*polyhedronB);
604  } else if (solid->GetEntityType() == "G4IntersectionSolid") {
605  *result = polyhedronA->intersect(*polyhedronB);
606  } else {
607  B2WARNING("getBooleanSolidPolyhedron(): Unrecognized boolean solid " << solid->GetName() <<
608  " of type " << solid->GetEntityType());
609  }
610  delete polyhedronA;
611  delete polyhedronB;
612  return result;
613 }
Belle2::Module::setDescription
void setDescription(const std::string &description)
Sets the description of the module.
Definition: Module.cc:214
Belle2::VRMLWriterModule::describeSolid
void describeSolid(G4VSolid *, const std::string &, bool)
Emit VRML for each solid.
Definition: VRMLWriterModule.cc:161
Belle2::VRMLWriterModule::VRMLWriterModule
VRMLWriterModule()
Constructor of the module.
Definition: VRMLWriterModule.cc:36
Belle2::VRMLWriterModule::writePolyhedron
void writePolyhedron(const G4Polyhedron *, const std::string &)
Emit VRML for the solid's polyhedron.
Definition: VRMLWriterModule.cc:483
Belle2::VRMLWriterModule::initialize
void initialize() override
Initialize at the start of a job.
Definition: VRMLWriterModule.cc:45
Belle2::VRMLWriterModule::event
void event() override
Called for each event: this runs the VRML writer only for the first event.
Definition: VRMLWriterModule.cc:50
Belle2::VRMLWriterModule::writePhysicalVolume
void writePhysicalVolume(const G4VPhysicalVolume *, const std::string &, const std::string &, bool)
Emit VRML for each daughter of a logical volume.
Definition: VRMLWriterModule.cc:516
Belle2::VRMLWriterModule::descendAndDescribe
void descendAndDescribe(G4VPhysicalVolume *, const std::string &, int)
Emit VRML for a physical volume (recursive)
Definition: VRMLWriterModule.cc:423
Belle2::VRMLWriterModule::m_File
std::ofstream m_File
Output file.
Definition: VRMLWriterModule.h:90
Belle2::VRMLWriterModule::writePreamble
void writePreamble(void)
Emit VRML for the start of the file.
Definition: VRMLWriterModule.cc:445
Belle2::VRMLWriterModule::m_First
bool m_First
Once-only flag to write VRML only on the first event.
Definition: VRMLWriterModule.h:84
Belle2::VRMLWriterModule::m_PVIndex
std::vector< std::vector< int > > * m_PVIndex
Indices (in G4PhysicalVolumeStore) of the logical volume's physical-volume daughters.
Definition: VRMLWriterModule.h:105
Belle2::VRMLWriterModule::assignName
void assignName(std::vector< std::string > *, unsigned int, const G4String &, int)
Create unique and legal name for each solid.
Definition: VRMLWriterModule.cc:126
Belle2::VRMLWriterModule::m_IsCylinder
std::vector< bool > * m_IsCylinder
Flag to indicate that a solid can be rendered as a VMRL cylinder.
Definition: VRMLWriterModule.h:102
Belle2::VRMLWriterModule::m_Filename
std::string m_Filename
User-specified output filename.
Definition: VRMLWriterModule.h:87
Belle2::VRMLWriterModule::m_PVWritten
std::vector< bool > * m_PVWritten
Flag to indicate that the physical volume has already been written.
Definition: VRMLWriterModule.h:111
Belle2::VRMLWriterModule::m_PVName
std::vector< std::string > * m_PVName
Modified (legal-character and unique) physical-volume name.
Definition: VRMLWriterModule.h:93
Belle2::VRMLWriterModule::m_SolidName
std::vector< std::string > * m_SolidName
Modified (legal-character and unique) solid name.
Definition: VRMLWriterModule.h:99
Belle2::VRMLWriterModule::describeLogicalVolume
void describeLogicalVolume(G4LogicalVolume *, const std::string &, const std::string &, bool)
Emit VRML for each logical volume.
Definition: VRMLWriterModule.cc:196
Belle2::VRMLWriterModule::m_LVWritten
std::vector< bool > * m_LVWritten
Flag to indicate that the logical volume has already been written.
Definition: VRMLWriterModule.h:108
Belle2::VRMLWriterModule::getBooleanSolidPolyhedron
HepPolyhedron * getBooleanSolidPolyhedron(G4VSolid *)
Create polyhedron for a boolean solid (recursive)
Definition: VRMLWriterModule.cc:568
Belle2::VRMLWriterModule::m_LVName
std::vector< std::string > * m_LVName
Modified (legal-character and unique) logical-volume name.
Definition: VRMLWriterModule.h:96
Belle2::VRMLWriterModule::describePhysicalVolume
void describePhysicalVolume(G4VPhysicalVolume *)
Access next physical volume in the tree (recursive)
Definition: VRMLWriterModule.cc:356
Belle2::geometry::GeometryManager::getInstance
static GeometryManager & getInstance()
Return a reference to the instance.
Definition: GeometryManager.cc:96
Belle2::geometry::GeometryManager::getTopVolume
G4VPhysicalVolume * getTopVolume()
Return a pointer to the top volume.
Definition: GeometryManager.h:56
Belle2::REG_MODULE
REG_MODULE(arichBtest)
Register the Module.
Belle2::Module::addParam
void addParam(const std::string &name, T &paramVariable, const std::string &description, const T &defaultValue)
Adds a new parameter to the module.
Definition: Module.h:560
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:17