Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
TGeoVGShape.cxx
Go to the documentation of this file.
1 // Author: Mihaela Gheata 30/03/16
2 /*************************************************************************
3  * Copyright (C) 1995-2016, Rene Brun and Fons Rademakers. *
4  * All rights reserved. *
5  * *
6  * For the licensing terms see $ROOTSYS/LICENSE. *
7  * For the list of contributors see $ROOTSYS/README/CREDITS. *
8  *************************************************************************/
9 
10 /** \class TGeoVGShape
11 \ingroup Geometry_classes
12 
13 Bridge class for using a VecGeom solid as TGeoShape.
14 */
15 
16 #include "TGeoVGShape.h"
17 
18 #include "volumes/PlacedVolume.h"
19 #include "volumes/UnplacedVolume.h"
20 #include "volumes/UnplacedBox.h"
21 #include "volumes/UnplacedTube.h"
22 #include "volumes/UnplacedCone.h"
23 #include "volumes/UnplacedParaboloid.h"
24 #include "volumes/UnplacedParallelepiped.h"
25 #include "volumes/UnplacedPolyhedron.h"
26 #include "volumes/UnplacedTrd.h"
27 #include "volumes/UnplacedOrb.h"
28 #include "volumes/UnplacedSphere.h"
29 #include "volumes/UnplacedBooleanVolume.h"
30 #include "volumes/UnplacedTorus2.h"
31 #include "volumes/UnplacedTrapezoid.h"
32 #include "volumes/UnplacedPolycone.h"
33 #include "volumes/UnplacedScaledShape.h"
34 #include "volumes/UnplacedGenTrap.h"
35 #include "volumes/UnplacedSExtruVolume.h"
36 #include "TError.h"
37 #include "TGeoManager.h"
38 #include "TGeoMaterial.h"
39 #include "TGeoMedium.h"
40 #include "TGeoVolume.h"
41 #include "TGeoArb8.h"
42 #include "TGeoTube.h"
43 #include "TGeoCone.h"
44 #include "TGeoPara.h"
45 #include "TGeoParaboloid.h"
46 #include "TGeoTrd1.h"
47 #include "TGeoTrd2.h"
48 #include "TGeoPcon.h"
49 #include "TGeoPgon.h"
50 #include "TGeoSphere.h"
51 #include "TGeoBoolNode.h"
52 #include "TGeoCompositeShape.h"
53 #include "TGeoScaledShape.h"
54 #include "TGeoTorus.h"
55 #include "TGeoEltu.h"
56 #include "TGeoXtru.h"
57 
58 ////////////////////////////////////////////////////////////////////////////////
59 /// Default constructor
60 
61 TGeoVGShape::TGeoVGShape(TGeoShape *shape, vecgeom::cxx::VPlacedVolume *vgshape)
62  :TGeoBBox(shape->GetName(), 0, 0, 0), fVGShape(vgshape), fShape(shape)
63 {
64  // Copy box parameters from the original ROOT shape
65  const TGeoBBox *box = (const TGeoBBox*)shape;
66  TGeoBBox::SetBoxDimensions(box->GetDX(), box->GetDY(), box->GetDZ());
67  memcpy(fOrigin, box->GetOrigin(), 3*sizeof(Double_t));
68 }
69 
70 ////////////////////////////////////////////////////////////////////////////////
71 /// Destructor
72 
73 TGeoVGShape::~TGeoVGShape()
74 {
75  // Cleanup only the VecGeom solid, the ROOT shape is cleaned by TGeoManager
76  delete fVGShape;
77 }
78 
79 ////////////////////////////////////////////////////////////////////////////////
80 /// Factory creating TGeoVGShape from a Root shape. Returns nullptr if the
81 /// shape cannot be converted
82 
83 TGeoVGShape *TGeoVGShape::Create(TGeoShape *shape)
84 {
85  vecgeom::cxx::VPlacedVolume *vgshape = TGeoVGShape::CreateVecGeomSolid(shape);
86  if (!vgshape) return nullptr;
87  return ( new TGeoVGShape(shape, vgshape) );
88 }
89 
90 ////////////////////////////////////////////////////////////////////////////////
91 /// Conversion method to create VecGeom solid corresponding to TGeoShape
92 
93 vecgeom::cxx::VPlacedVolume *TGeoVGShape::CreateVecGeomSolid(TGeoShape *shape)
94 {
95  // Call VecGeom TGeoShape->UnplacedSolid converter
96  // VUnplacedVolume *unplaced = RootGeoManager::Instance().Convert(shape);
97  vecgeom::cxx::VUnplacedVolume *unplaced = Convert(shape);
98  if (!unplaced) return nullptr;
99  // We have to create a placed volume from the unplaced one to have access
100  // to the navigation interface
101  vecgeom::cxx::LogicalVolume *lvol = new vecgeom::cxx::LogicalVolume("", unplaced);
102  return ( lvol->Place() );
103 }
104 
105 ////////////////////////////////////////////////////////////////////////////////
106 /// Convert a TGeoMatrix to a TRansformation3D
107 
108 vecgeom::cxx::Transformation3D *TGeoVGShape::Convert(TGeoMatrix const *const geomatrix) {
109  Double_t const *const t = geomatrix->GetTranslation();
110  Double_t const *const r = geomatrix->GetRotationMatrix();
111  vecgeom::cxx::Transformation3D *const transformation =
112  new vecgeom::cxx::Transformation3D(t[0], t[1], t[2], r[0], r[1], r[2], r[3], r[4], r[5], r[6], r[7], r[8]);
113  return transformation;
114 }
115 
116 ////////////////////////////////////////////////////////////////////////////////
117 /// Convert a TGeo shape to VUnplacedVolume, then creates a VPlacedVolume
118 
119 vecgeom::cxx::VUnplacedVolume* TGeoVGShape::Convert(TGeoShape const *const shape)
120 {
121  using namespace vecgeom::cxx;
122  VUnplacedVolume *unplaced_volume = nullptr;
123 
124  // THE BOX
125  if (shape->IsA() == TGeoBBox::Class()) {
126  TGeoBBox const *const box = static_cast<TGeoBBox const *>(shape);
127  unplaced_volume = new UnplacedBox(box->GetDX(), box->GetDY(), box->GetDZ());
128  }
129 
130  // THE TUBE
131  if (shape->IsA() == TGeoTube::Class()) {
132  TGeoTube const *const tube = static_cast<TGeoTube const *>(shape);
133  unplaced_volume = new GenericUnplacedTube(tube->GetRmin(), tube->GetRmax(), tube->GetDz(), 0., kTwoPi);
134  }
135 
136  // THE TUBESEG
137  if (shape->IsA() == TGeoTubeSeg::Class()) {
138  TGeoTubeSeg const *const tube = static_cast<TGeoTubeSeg const *>(shape);
139  unplaced_volume =
140  new GenericUnplacedTube(tube->GetRmin(), tube->GetRmax(), tube->GetDz(), kDegToRad * tube->GetPhi1(),
141  kDegToRad * (tube->GetPhi2() - tube->GetPhi1()));
142  }
143 
144  // THE CONESEG
145  if (shape->IsA() == TGeoConeSeg::Class()) {
146  TGeoConeSeg const *const cone = static_cast<TGeoConeSeg const *>(shape);
147  unplaced_volume =
148  new UnplacedCone(cone->GetRmin1(), cone->GetRmax1(), cone->GetRmin2(), cone->GetRmax2(), cone->GetDz(),
149  kDegToRad * cone->GetPhi1(), kDegToRad * (cone->GetPhi2() - cone->GetPhi1()));
150  }
151 
152  // THE CONE
153  if (shape->IsA() == TGeoCone::Class()) {
154  TGeoCone const *const cone = static_cast<TGeoCone const *>(shape);
155  unplaced_volume = new UnplacedCone(cone->GetRmin1(), cone->GetRmax1(), cone->GetRmin2(), cone->GetRmax2(),
156  cone->GetDz(), 0., kTwoPi);
157  }
158 
159  // THE PARABOLOID
160  if (shape->IsA() == TGeoParaboloid::Class()) {
161  TGeoParaboloid const *const p = static_cast<TGeoParaboloid const *>(shape);
162  unplaced_volume = new UnplacedParaboloid(p->GetRlo(), p->GetRhi(), p->GetDz());
163  }
164 
165  // THE PARALLELEPIPED
166  if (shape->IsA() == TGeoPara::Class()) {
167  TGeoPara const *const p = static_cast<TGeoPara const *>(shape);
168  unplaced_volume =
169  new UnplacedParallelepiped(p->GetX(), p->GetY(), p->GetZ(), p->GetAlpha(), p->GetTheta(), p->GetPhi());
170  }
171 
172  // Polyhedron/TGeoPgon
173  if (shape->IsA() == TGeoPgon::Class()) {
174  TGeoPgon const *pgon = static_cast<TGeoPgon const *>(shape);
175  unplaced_volume = new UnplacedPolyhedron(pgon->GetPhi1(), // phiStart
176  pgon->GetDphi(), // phiEnd
177  pgon->GetNedges(), // sideCount
178  pgon->GetNz(), // zPlaneCount
179  pgon->GetZ(), // zPlanes
180  pgon->GetRmin(), // rMin
181  pgon->GetRmax() // rMax
182  );
183  }
184 
185  // TRD2
186  if (shape->IsA() == TGeoTrd2::Class()) {
187  TGeoTrd2 const *const p = static_cast<TGeoTrd2 const *>(shape);
188  unplaced_volume = new UnplacedTrd(p->GetDx1(), p->GetDx2(), p->GetDy1(), p->GetDy2(), p->GetDz());
189  }
190 
191  // TRD1
192  if (shape->IsA() == TGeoTrd1::Class()) {
193  TGeoTrd1 const *const p = static_cast<TGeoTrd1 const *>(shape);
194  unplaced_volume = new UnplacedTrd(p->GetDx1(), p->GetDx2(), p->GetDy(), p->GetDz());
195  }
196 
197  // TRAPEZOID
198  if (shape->IsA() == TGeoTrap::Class()) {
199  TGeoTrap const *const p = static_cast<TGeoTrap const *>(shape);
200  unplaced_volume = new UnplacedTrapezoid(p->GetDz(), p->GetTheta() * kDegToRad, p->GetPhi() * kDegToRad, p->GetH1(),
201  p->GetBl1(), p->GetTl1(), std::tan(p->GetAlpha1() * kDegToRad), p->GetH2(),
202  p->GetBl2(), p->GetTl2(), std::tan(p->GetAlpha2() * kDegToRad));
203  }
204 
205  // THE SPHERE | ORB
206  if (shape->IsA() == TGeoSphere::Class()) {
207  // make distinction
208  TGeoSphere const *const p = static_cast<TGeoSphere const *>(shape);
209  if (p->GetRmin() == 0. && p->GetTheta2() - p->GetTheta1() == 180. && p->GetPhi2() - p->GetPhi1() == 360.) {
210  unplaced_volume = new UnplacedOrb(p->GetRmax());
211  } else {
212  unplaced_volume = new UnplacedSphere(p->GetRmin(), p->GetRmax(), p->GetPhi1() * kDegToRad,
213  (p->GetPhi2() - p->GetPhi1()) * kDegToRad, p->GetTheta1() * kDegToRad,
214  (p->GetTheta2() - p->GetTheta1()) * kDegToRad);
215  }
216  }
217 
218  if (shape->IsA() == TGeoCompositeShape::Class()) {
219  TGeoCompositeShape const *const compshape = static_cast<TGeoCompositeShape const *>(shape);
220  TGeoBoolNode const *const boolnode = compshape->GetBoolNode();
221 
222  // need the matrix;
223  Transformation3D const *lefttrans = Convert(boolnode->GetLeftMatrix());
224  Transformation3D const *righttrans = Convert(boolnode->GetRightMatrix());
225  // unplaced shapes
226  VUnplacedVolume const *leftunplaced = Convert(boolnode->GetLeftShape());
227  VUnplacedVolume const *rightunplaced = Convert(boolnode->GetRightShape());
228  if (!leftunplaced || !rightunplaced) {
229  // If one of the components cannot be converted, cleanup & return nullptr
230  delete lefttrans;
231  delete righttrans;
232  delete leftunplaced;
233  delete rightunplaced;
234  return nullptr;
235  }
236 
237  assert(leftunplaced != nullptr);
238  assert(rightunplaced != nullptr);
239 
240  // the problem is that I can only place logical volumes
241  VPlacedVolume *const leftplaced = (new LogicalVolume("inner_virtual", leftunplaced))->Place(lefttrans);
242 
243  VPlacedVolume *const rightplaced = (new LogicalVolume("inner_virtual", rightunplaced))->Place(righttrans);
244 
245  // now it depends on concrete type
246  if (boolnode->GetBooleanOperator() == TGeoBoolNode::kGeoSubtraction) {
247  unplaced_volume = new UnplacedBooleanVolume(kSubtraction, leftplaced, rightplaced);
248  } else if (boolnode->GetBooleanOperator() == TGeoBoolNode::kGeoIntersection) {
249  unplaced_volume = new UnplacedBooleanVolume(kIntersection, leftplaced, rightplaced);
250  } else if (boolnode->GetBooleanOperator() == TGeoBoolNode::kGeoUnion) {
251  unplaced_volume = new UnplacedBooleanVolume(kUnion, leftplaced, rightplaced);
252  }
253  }
254 
255  // THE TORUS
256  if (shape->IsA() == TGeoTorus::Class()) {
257  // make distinction
258  TGeoTorus const *const p = static_cast<TGeoTorus const *>(shape);
259  unplaced_volume =
260  new UnplacedTorus2(p->GetRmin(), p->GetRmax(), p->GetR(), p->GetPhi1() * kDegToRad, p->GetDphi() * kDegToRad);
261  }
262 
263  // THE POLYCONE
264  if (shape->IsA() == TGeoPcon::Class()) {
265  TGeoPcon const *const p = static_cast<TGeoPcon const *>(shape);
266  unplaced_volume = new UnplacedPolycone(p->GetPhi1() * kDegToRad, p->GetDphi() * kDegToRad, p->GetNz(), p->GetZ(),
267  p->GetRmin(), p->GetRmax());
268  }
269 
270  // THE SCALED SHAPE
271  if (shape->IsA() == TGeoScaledShape::Class()) {
272  TGeoScaledShape const *const p = static_cast<TGeoScaledShape const *>(shape);
273  // First convert the referenced shape
274  VUnplacedVolume *referenced_shape = Convert(p->GetShape());
275  if (!referenced_shape) return nullptr;
276  const double *scale_root = p->GetScale()->GetScale();
277  unplaced_volume = new UnplacedScaledShape(referenced_shape, scale_root[0], scale_root[1], scale_root[2]);
278  }
279 
280  // THE ELLIPTICAL TUBE AS SCALED TUBE
281  if (shape->IsA() == TGeoEltu::Class()) {
282  TGeoEltu const *const p = static_cast<TGeoEltu const *>(shape);
283  // Create the corresponding unplaced tube, with:
284  // rmin=0, rmax=A, dz=dz, which is scaled with (1., A/B, 1.)
285  GenericUnplacedTube *tubeUnplaced = new GenericUnplacedTube(0, p->GetA(), p->GetDZ(), 0, kTwoPi);
286  unplaced_volume = new UnplacedScaledShape(tubeUnplaced, 1., p->GetB() / p->GetA(), 1.);
287  }
288 
289  // THE ARB8
290  if (shape->IsA() == TGeoArb8::Class() || shape->IsA() == TGeoGtra::Class()) {
291  TGeoArb8 *p = (TGeoArb8 *)(shape);
292  // Create the corresponding GenTrap
293  std::vector<Vector3D<Precision>> vertexlist;
294  const double *vertices = p->GetVertices();
295  Precision verticesx[8], verticesy[8];
296  for (auto ivert = 0; ivert < 8; ++ivert) {
297  verticesx[ivert] = vertices[2 * ivert];
298  verticesy[ivert] = vertices[2 * ivert + 1];
299  }
300  unplaced_volume = new UnplacedGenTrap(verticesx, verticesy, p->GetDz());
301  }
302 
303  // THE SIMPLE XTRU
304  if (shape->IsA() == TGeoXtru::Class()) {
305  TGeoXtru *p = (TGeoXtru *)(shape);
306  // analyse convertibility
307  if (p->GetNz() == 2) {
308  // add check on scaling and distortions
309  size_t Nvert = (size_t)p->GetNvert();
310  double *x = new double[Nvert];
311  double *y = new double[Nvert];
312  for (size_t i = 0; i < Nvert; ++i) {
313  x[i] = p->GetX(i);
314  y[i] = p->GetY(i);
315  }
316  // check in which orientation the polygon in given
317  if (PlanarPolygon::GetOrientation(x, y, Nvert) > 0.) {
318  // std::cerr << "Points not given in clockwise order ... reordering \n";
319  for (size_t i = 0; i < Nvert; ++i) {
320  x[Nvert - 1 - i] = p->GetX(i);
321  y[Nvert - 1 - i] = p->GetY(i);
322  }
323  }
324  unplaced_volume = new UnplacedSExtruVolume(p->GetNvert(), x, y, p->GetZ()[0], p->GetZ()[1]);
325  }
326  }
327 
328  // New volumes should be implemented here...
329  if (!unplaced_volume) {
330  printf("Unsupported shape for ROOT shape \"%s\" of type %s. "
331  "Using ROOT implementation.\n",
332  shape->GetName(), shape->ClassName());
333  return nullptr;
334  }
335 
336  return ( unplaced_volume );
337 }
338 
339 ////////////////////////////////////////////////////////////////////////////////
340 /// Compute bounding box.
341 
342 void TGeoVGShape::ComputeBBox()
343 {
344  fShape->ComputeBBox();
345 }
346 
347 ////////////////////////////////////////////////////////////////////////////////
348 /// Returns analytic capacity of the solid
349 
350 Double_t TGeoVGShape::Capacity() const
351 {
352  return fVGShape->Capacity();
353 }
354 
355 ////////////////////////////////////////////////////////////////////////////////
356 /// Normal computation.
357 
358 void TGeoVGShape::ComputeNormal(const Double_t *point, const Double_t */*dir*/, Double_t *norm)
359 {
360  vecgeom::cxx::Vector3D<Double_t> vnorm;
361  fVGShape->Normal(vecgeom::cxx::Vector3D<Double_t>(point[0], point[1], point[2]), vnorm);
362  norm[0] = vnorm.x(); norm[1] = vnorm.y(), norm[2] = vnorm.z();
363 }
364 
365 ////////////////////////////////////////////////////////////////////////////////
366 /// Test if point is inside this shape.
367 
368 Bool_t TGeoVGShape::Contains(const Double_t *point) const
369 {
370  return ( fVGShape->Contains(vecgeom::cxx::Vector3D<Double_t>(point[0], point[1], point[2])) );
371 }
372 
373 ////////////////////////////////////////////////////////////////////////////////
374 
375 Double_t TGeoVGShape::DistFromInside(const Double_t *point, const Double_t *dir, Int_t /*iact*/,
376  Double_t step, Double_t * /*safe*/) const
377 {
378  Double_t dist = fVGShape->DistanceToOut(vecgeom::cxx::Vector3D<Double_t>(point[0], point[1], point[2]),
379  vecgeom::cxx::Vector3D<Double_t>(dir[0], dir[1], dir[2]), step);
380  return ( (dist < 0.)? 0. : dist );
381 }
382 
383 ////////////////////////////////////////////////////////////////////////////////
384 
385 Double_t TGeoVGShape::DistFromOutside(const Double_t *point, const Double_t *dir, Int_t /*iact*/,
386  Double_t step, Double_t * /*safe*/) const
387 {
388  Double_t dist = fVGShape->DistanceToIn(vecgeom::cxx::Vector3D<Double_t>(point[0], point[1], point[2]),
389  vecgeom::cxx::Vector3D<Double_t>(dir[0], dir[1], dir[2]), step);
390  return ( (dist < 0.)? 0. : dist );
391 }
392 
393 ////////////////////////////////////////////////////////////////////////////////
394 
395 Double_t TGeoVGShape::Safety(const Double_t *point, Bool_t in) const
396 {
397  Double_t safety = (in) ? fVGShape->SafetyToOut(vecgeom::cxx::Vector3D<Double_t>(point[0], point[1], point[2]))
398  : fVGShape->SafetyToIn(vecgeom::cxx::Vector3D<Double_t>(point[0], point[1], point[2]));
399  return ( (safety < 0.)? 0. : safety );
400 }
401 
402 ////////////////////////////////////////////////////////////////////////////////
403 /// Print info about the VecGeom solid
404 
405 void TGeoVGShape::InspectShape() const
406 {
407  fVGShape->GetUnplacedVolume()->Print();
408  printf("\n");
409 }