Logo ROOT  

Reference Guide
Loading...
Searching...
No Matches
TGeoArb8.cxx
Go to the documentation of this file.
1// @(#)root/geom:$Id$
2// Author: Andrei Gheata 31/01/02
3
4/*************************************************************************
5 * Copyright (C) 1995-2000, 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#include "TGeoArb8.h"
13
14#include <iostream>
15#include "TBuffer.h"
16#include "TGeoManager.h"
17#include "TGeoVolume.h"
18#include "TGeoMatrix.h"
19#include "TMath.h"
20
21
22/** \class TGeoArb8
23\ingroup Trapezoids
24
25\brief An arbitrary trapezoid with less than 8 vertices standing on two parallel planes perpendicular to Z axis.
26
27
28An `Arb8` is defined by two quadrilaterals sitting on parallel planes,
29at `dZ`. These are defined each by 4 vertices having the coordinates
30`(Xi,Yi,+/-dZ)`,` i=0`,` 3`. The lateral surface of the `Arb8` is
31defined by the 4 pairs of edges corresponding to vertices (`i,i+1`) on
32both `-dZ` and `+dZ`. If M and M' are the middles of the segments
33`(i,i+1)` at `-dZ` and `+dZ`, a lateral surface is obtained by sweeping
34the edge at `-dZ` along MM' so that it will match the corresponding one
35at `+dZ`. Since the points defining the edges are arbitrary, the lateral
36surfaces are not necessary planes - but twisted planes having a twist
37angle linear-dependent on Z.
38
39~~~ {.cpp}
40TGeoArb8::TGeoArb8(Double_t dz,Double_t ivert);
41~~~
42
43 - `dz:` half-length in Z;
44 - `ivert = [0,7]`
45
46Vertices have to be defined clockwise in the XY pane, both at `+dz` and
47`-dz`. The quadrilateral at `-dz` is defined by indices [0,3], whereas
48the one at `+dz` by vertices [4,7]. The vertex with `index=7` has to be
49defined last, since it triggers the computation of the bounding box of
50the shape. Any two or more vertices in each Z plane can have the same
51(X,Y) coordinates. It this case, the top and bottom quadrilaterals
52become triangles, segments or points. The lateral surfaces are not
53necessary defined by a pair of segments, but by pair segment-point
54(making a triangle) or point-point (making a line). Any choice is valid
55as long as at one of the end-caps is at least a triangle.
56
57Begin_Macro
58{
59 TCanvas *c = new TCanvas("c", "c",0,0,600,600);
60 new TGeoManager("arb8", "poza12");
61 TGeoMaterial *mat = new TGeoMaterial("Al", 26.98,13,2.7);
62 TGeoMedium *med = new TGeoMedium("MED",1,mat);
63 TGeoVolume *top = gGeoManager->MakeBox("TOP",med,100,100,100);
64 gGeoManager->SetTopVolume(top);
65 TGeoArb8 *arb = new TGeoArb8(20);
66 arb->SetVertex(0,-30,-25);
67 arb->SetVertex(1,-25,25);
68 arb->SetVertex(2,5,25);
69 arb->SetVertex(3,25,-25);
70 arb->SetVertex(4,-28,-23);
71 arb->SetVertex(5,-23,27);
72 arb->SetVertex(6,-23,27);
73 arb->SetVertex(7,13,-27);
74 TGeoVolume *vol = new TGeoVolume("ARB8",arb,med);
75 vol->SetLineWidth(2);
76 top->AddNode(vol,1);
77 gGeoManager->CloseGeometry();
78 gGeoManager->SetNsegments(80);
79 top->Draw();
80 TView *view = gPad->GetView();
81 if (view) view->ShowAxis();
82}
83End_Macro
84*/
85
86/** \class TGeoGtra
87\ingroup Trapezoids
88\brief A twisted trapezoid.
89
90A twisted trapezoid is a general trapezoid defined in the same way but
91that is twisted along the Z-axis. The twist is defined as the rotation
92angle between the lower and the higher Z faces.
93
94~~~ {.cpp}
95TGeoGtra(Double_t dz,Double_t theta,Double_t phi,Double_t twist,
96Double_t h1,Double_t bl1,Double_t tl1,Double_t alpha1,
97Double_t h2,Double_t bl2,Double_t tl2,Double_t alpha2 );
98~~~
99
100Begin_Macro
101{
102 TCanvas *c = new TCanvas("c", "c",0,0,600,600);
103 new TGeoManager("gtra", "poza11");
104 TGeoMaterial *mat = new TGeoMaterial("Al", 26.98,13,2.7);
105 TGeoMedium *med = new TGeoMedium("MED",1,mat);
106 TGeoVolume *top = gGeoManager->MakeBox("TOP",med,100,100,100);
107 gGeoManager->SetTopVolume(top);
108 TGeoVolume *vol = gGeoManager->MakeGtra("Gtra",med, 30,15,30,30,20,10,15,0,20,10,15,0);
109 vol->SetLineWidth(2);
110 top->AddNode(vol,1);
111 gGeoManager->CloseGeometry();
112 gGeoManager->SetNsegments(80);
113 top->Draw();
114 TView *view = gPad->GetView();
115 if (view) view->ShowAxis();
116}
117End_Macro
118*/
119
120/** \class TGeoTrap
121\ingroup Trapezoids
122\brief A general trapezoid.
123
124A general trapezoid is one for which the faces perpendicular to z are
125trapezes but their centers are not necessary at the same x, y
126coordinates.
127
128Begin_Macro
129{
130 TCanvas *c = new TCanvas("c", "c",0,0,600,600);
131 new TGeoManager("trap", "poza10");
132 TGeoMaterial *mat = new TGeoMaterial("Al", 26.98,13,2.7);
133 TGeoMedium *med = new TGeoMedium("MED",1,mat);
134 TGeoVolume *top = gGeoManager->MakeBox("TOP",med,100,100,100);
135 gGeoManager->SetTopVolume(top);
136 TGeoVolume *vol = gGeoManager->MakeTrap("Trap",med, 30,15,30,20,10,15,0,20,10,15,0);
137 vol->SetLineWidth(2);
138 top->AddNode(vol,1);
139 gGeoManager->CloseGeometry();
140 gGeoManager->SetNsegments(80);
141 top->Draw();
142 TView *view = gPad->GetView();
143 if (view) view->ShowAxis();
144}
145End_Macro
146
147It has eleven parameters: the half length in z, the polar angles from
148the center of the face at low z to that at high z, `H1` the half length
149in y at low z, `LB1` the half length in x at low z and y low edge, `LB2`
150the half length in x at low z and y high edge, **`TH1`** the angle with
151respect to the y axis from the center of low y edge to the center of the
152high y edge, and `H2,LB2,LH2,TH2`, the corresponding quantities at high
153z.
154
155~~~ {.cpp}
156TGeoTrap(Double_t dz,Double_t theta,Double_t phi,
157Double_t h1,Double_t bl1,Double_t tl1,Double_t alpha1,
158Double_t h2,Double_t bl2,Double_t tl2,Double_t alpha2);
159~~~
160*/
161
162////////////////////////////////////////////////////////////////////////////////
163/// Default constructor.
164
165TGeoArb8::TGeoArb8()
166{
167 fDz = 0;
168 for (Int_t i = 0; i < 8; i++) {
169 fXY[i][0] = 0.0;
170 fXY[i][1] = 0.0;
171 }
172 SetShapeBit(kGeoArb8);
173}
174
175////////////////////////////////////////////////////////////////////////////////
176/// Constructor. If the array of vertices is not null, this should be
177/// in the format : (x0, y0, x1, y1, ... , x7, y7)
178
179TGeoArb8::TGeoArb8(Double_t dz, Double_t *vertices) : TGeoBBox(0, 0, 0)
180{
181 fDz = dz;
182 SetShapeBit(kGeoArb8);
183 if (vertices) {
184 for (Int_t i = 0; i < 8; i++) {
185 fXY[i][0] = vertices[2 * i];
186 fXY[i][1] = vertices[2 * i + 1];
187 }
188 ComputeTwist();
189 ComputeBBox();
190 } else {
191 for (Int_t i = 0; i < 8; i++) {
192 fXY[i][0] = 0.0;
193 fXY[i][1] = 0.0;
194 }
195 }
196}
197
198////////////////////////////////////////////////////////////////////////////////
199/// Named constructor. If the array of vertices is not null, this should be
200/// in the format : (x0, y0, x1, y1, ... , x7, y7)
201
202TGeoArb8::TGeoArb8(const char *name, Double_t dz, Double_t *vertices) : TGeoBBox(name, 0, 0, 0)
203{
204 fDz = dz;
205 SetShapeBit(kGeoArb8);
206 if (vertices) {
207 for (Int_t i = 0; i < 8; i++) {
208 fXY[i][0] = vertices[2 * i];
209 fXY[i][1] = vertices[2 * i + 1];
210 }
211 ComputeTwist();
212 ComputeBBox();
213 } else {
214 for (Int_t i = 0; i < 8; i++) {
215 fXY[i][0] = 0.0;
216 fXY[i][1] = 0.0;
217 }
218 }
219}
220
221////////////////////////////////////////////////////////////////////////////////
222/// Destructor.
223
224TGeoArb8::~TGeoArb8()
225{
226 if (fTwist)
227 delete[] fTwist;
228}
229
230////////////////////////////////////////////////////////////////////////////////
231/// Copy twist values from source array
232
233void TGeoArb8::CopyTwist(Double_t *twist)
234{
235 if (twist) {
236 if (!fTwist)
237 fTwist = new Double_t[4];
238 memcpy(fTwist, twist, 4 * sizeof(Double_t));
239 } else if (fTwist) {
240 delete[] fTwist;
241 fTwist = nullptr;
242 }
243}
244
245////////////////////////////////////////////////////////////////////////////////
246/// Computes capacity of the shape in [length^3].
247
248Double_t TGeoArb8::Capacity() const
249{
250 Int_t i, j;
251 Double_t capacity = 0;
252 for (i = 0; i < 4; i++) {
253 j = (i + 1) % 4;
254 capacity += 0.25 * fDz *
255 ((fXY[i][0] + fXY[i + 4][0]) * (fXY[j][1] + fXY[j + 4][1]) -
256 (fXY[j][0] + fXY[j + 4][0]) * (fXY[i][1] + fXY[i + 4][1]) +
257 (1. / 3) * ((fXY[i + 4][0] - fXY[i][0]) * (fXY[j + 4][1] - fXY[j][1]) -
258 (fXY[j][0] - fXY[j + 4][0]) * (fXY[i][1] - fXY[i + 4][1])));
259 }
260 return TMath::Abs(capacity);
261}
262
263////////////////////////////////////////////////////////////////////////////////
264/// Computes bounding box for an Arb8 shape.
265
266void TGeoArb8::ComputeBBox()
267{
268 Double_t xmin, xmax, ymin, ymax;
269 xmin = xmax = fXY[0][0];
270 ymin = ymax = fXY[0][1];
271
272 for (Int_t i = 1; i < 8; i++) {
273 if (xmin > fXY[i][0])
274 xmin = fXY[i][0];
275 if (xmax < fXY[i][0])
276 xmax = fXY[i][0];
277 if (ymin > fXY[i][1])
278 ymin = fXY[i][1];
279 if (ymax < fXY[i][1])
280 ymax = fXY[i][1];
281 }
282 fDX = 0.5 * (xmax - xmin);
283 fDY = 0.5 * (ymax - ymin);
284 fDZ = fDz;
285 fOrigin[0] = 0.5 * (xmax + xmin);
286 fOrigin[1] = 0.5 * (ymax + ymin);
287 fOrigin[2] = 0;
288 SetShapeBit(kGeoClosedShape);
289}
290
291////////////////////////////////////////////////////////////////////////////////
292/// Computes tangents of twist angles (angles between projections on XY plane
293/// of corresponding -dz +dz edges). Computes also if the vertices are defined
294/// clockwise or anti-clockwise.
295
296void TGeoArb8::ComputeTwist()
297{
298 Double_t twist[4];
299 Bool_t twisted = kFALSE;
300 Double_t dx1, dy1, dx2, dy2;
301 Bool_t singleBottom = kTRUE;
302 Bool_t singleTop = kTRUE;
303 Int_t i;
304 for (i = 0; i < 4; i++) {
305 dx1 = fXY[(i + 1) % 4][0] - fXY[i][0];
306 dy1 = fXY[(i + 1) % 4][1] - fXY[i][1];
307 if (TMath::Abs(dx1) < TGeoShape::Tolerance() && TMath::Abs(dy1) < TGeoShape::Tolerance()) {
308 twist[i] = 0;
309 continue;
310 }
311 singleBottom = kFALSE;
312 dx2 = fXY[4 + (i + 1) % 4][0] - fXY[4 + i][0];
313 dy2 = fXY[4 + (i + 1) % 4][1] - fXY[4 + i][1];
314 if (TMath::Abs(dx2) < TGeoShape::Tolerance() && TMath::Abs(dy2) < TGeoShape::Tolerance()) {
315 twist[i] = 0;
316 continue;
317 }
318 singleTop = kFALSE;
319 twist[i] = dy1 * dx2 - dx1 * dy2;
320 if (TMath::Abs(twist[i]) < TGeoShape::Tolerance()) {
321 twist[i] = 0;
322 continue;
323 }
324 twist[i] = TMath::Sign(1., twist[i]);
325 twisted = kTRUE;
326 }
327
328 CopyTwist(twisted ? twist : nullptr);
329
330 if (singleBottom) {
331 for (i = 0; i < 4; i++) {
332 fXY[i][0] += 1.E-8 * fXY[i + 4][0];
333 fXY[i][1] += 1.E-8 * fXY[i + 4][1];
334 }
335 }
336 if (singleTop) {
337 for (i = 0; i < 4; i++) {
338 fXY[i + 4][0] += 1.E-8 * fXY[i][0];
339 fXY[i + 4][1] += 1.E-8 * fXY[i][1];
340 }
341 }
342 Double_t sum1 = 0.;
343 Double_t sum2 = 0.;
344 Int_t j;
345 for (i = 0; i < 4; i++) {
346 j = (i + 1) % 4;
347 sum1 += fXY[i][0] * fXY[j][1] - fXY[j][0] * fXY[i][1];
348 sum2 += fXY[i + 4][0] * fXY[j + 4][1] - fXY[j + 4][0] * fXY[i + 4][1];
349 }
350 if (sum1 * sum2 < -TGeoShape::Tolerance()) {
351 Fatal("ComputeTwist", "Shape %s type Arb8: Lower/upper faces defined with opposite clockwise", GetName());
352 return;
353 }
354 if (sum1 > TGeoShape::Tolerance()) {
355 Error("ComputeTwist", "Shape %s type Arb8: Vertices must be defined clockwise in XY planes. Re-ordering...",
356 GetName());
357 Double_t xtemp, ytemp;
358 xtemp = fXY[1][0];
359 ytemp = fXY[1][1];
360 fXY[1][0] = fXY[3][0];
361 fXY[1][1] = fXY[3][1];
362 fXY[3][0] = xtemp;
363 fXY[3][1] = ytemp;
364 xtemp = fXY[5][0];
365 ytemp = fXY[5][1];
366 fXY[5][0] = fXY[7][0];
367 fXY[5][1] = fXY[7][1];
368 fXY[7][0] = xtemp;
369 fXY[7][1] = ytemp;
370 }
371 // Check for illegal crossings.
372 Bool_t illegal_cross = kFALSE;
373 illegal_cross =
374 TGeoShape::IsSegCrossing(fXY[0][0], fXY[0][1], fXY[1][0], fXY[1][1], fXY[2][0], fXY[2][1], fXY[3][0], fXY[3][1]);
375 if (!illegal_cross)
376 illegal_cross = TGeoShape::IsSegCrossing(fXY[4][0], fXY[4][1], fXY[5][0], fXY[5][1], fXY[6][0], fXY[6][1],
377 fXY[7][0], fXY[7][1]);
378 if (illegal_cross) {
379 Error("ComputeTwist", "Shape %s type Arb8: Malformed polygon with crossing opposite segments", GetName());
380 InspectShape();
381 }
382}
383
384////////////////////////////////////////////////////////////////////////////////
385/// Get twist for segment I in range [0,3]
386
387Double_t TGeoArb8::GetTwist(Int_t iseg) const
388{
389 return (!fTwist || iseg < 0 || iseg > 3) ? 0. : fTwist[iseg];
390}
391
392////////////////////////////////////////////////////////////////////////////////
393/// Get index of the edge of the quadrilater represented by vert closest to point.
394/// If [P1,P2] is the closest segment and P is the point, the function returns the fraction of the
395/// projection of (P1P) over (P1P2). If projection of P is not in range [P1,P2] return -1.
396
397Double_t TGeoArb8::GetClosestEdge(const Double_t *point, Double_t *vert, Int_t &isegment) const
398{
399 isegment = 0;
400 Int_t isegmin = 0;
401 Int_t i1, i2;
402 Double_t p1[2], p2[2];
403 Double_t lsq, ssq, dx, dy, dpx, dpy, u;
404 Double_t umin = -1.;
405 Double_t safe = 1E30;
406 for (i1 = 0; i1 < 4; i1++) {
407 if (TGeoShape::IsSameWithinTolerance(safe, 0)) {
408 isegment = isegmin;
409 return umin;
410 }
411 i2 = (i1 + 1) % 4;
412 p1[0] = vert[2 * i1];
413 p1[1] = vert[2 * i1 + 1];
414 p2[0] = vert[2 * i2];
415 p2[1] = vert[2 * i2 + 1];
416 dx = p2[0] - p1[0];
417 dy = p2[1] - p1[1];
418 dpx = point[0] - p1[0];
419 dpy = point[1] - p1[1];
420 lsq = dx * dx + dy * dy;
421 // Current segment collapsed to a point
422 if (TGeoShape::IsSameWithinTolerance(lsq, 0)) {
423 ssq = dpx * dpx + dpy * dpy;
424 if (ssq < safe) {
425 safe = ssq;
426 isegmin = i1;
427 umin = -1;
428 }
429 continue;
430 }
431 // Projection fraction
432 u = (dpx * dx + dpy * dy) / lsq;
433 // If projection of P is outside P1P2 limits, take the distance from P to the closest of P1 and P2
434 if (u > 1) {
435 // Outside, closer to P2
436 dpx = point[0] - p2[0];
437 dpy = point[1] - p2[1];
438 u = -1.;
439 } else {
440 if (u >= 0) {
441 // Projection inside
442 dpx -= u * dx;
443 dpy -= u * dy;
444 } else {
445 // Outside, closer to P1
446 u = -1.;
447 }
448 }
449 ssq = dpx * dpx + dpy * dpy;
450 if (ssq < safe) {
451 safe = ssq;
452 isegmin = i1;
453 umin = u;
454 }
455 }
456 isegment = isegmin;
457 // safe = TMath::Sqrt(safe);
458 return umin;
459}
460
461////////////////////////////////////////////////////////////////////////////////
462/// Compute normal to closest surface from POINT.
463
464void TGeoArb8::ComputeNormal(const Double_t *point, const Double_t *dir, Double_t *norm) const
465{
466 Double_t safc;
467 Double_t x0, y0, z0, x1, y1, z1, x2, y2, z2;
468 Double_t ax, ay, az, bx, by, bz;
469 Double_t fn;
470 safc = fDz - TMath::Abs(point[2]);
471 if (safc < 10. * TGeoShape::Tolerance()) {
472 memset(norm, 0, 3 * sizeof(Double_t));
473 norm[2] = (dir[2] > 0) ? 1 : (-1);
474 return;
475 }
476 Double_t vert[8];
477 SetPlaneVertices(point[2], vert);
478 // Get the closest edge (point should be on this edge within tolerance)
479 Int_t iseg;
480 Double_t frac = GetClosestEdge(point, vert, iseg);
481 if (frac < 0)
482 frac = 0.;
483 Int_t jseg = (iseg + 1) % 4;
484 x0 = vert[2 * iseg];
485 y0 = vert[2 * iseg + 1];
486 z0 = point[2];
487 x2 = vert[2 * jseg];
488 y2 = vert[2 * jseg + 1];
489 z2 = point[2];
490 x0 += frac * (x2 - x0);
491 y0 += frac * (y2 - y0);
492 x1 = fXY[iseg + 4][0];
493 y1 = fXY[iseg + 4][1];
494 z1 = fDz;
495 x1 += frac * (fXY[jseg + 4][0] - x1);
496 y1 += frac * (fXY[jseg + 4][1] - y1);
497 ax = x1 - x0;
498 ay = y1 - y0;
499 az = z1 - z0;
500 bx = x2 - x0;
501 by = y2 - y0;
502 bz = z2 - z0;
503 // Cross product of the vector given by the section segment (that contains the point) at z=point[2]
504 // and the vector connecting the point projection to its correspondent on the top edge (hard to swallow, isn't it?)
505 norm[0] = ay * bz - az * by;
506 norm[1] = az * bx - ax * bz;
507 norm[2] = ax * by - ay * bx;
508 fn = TMath::Sqrt(norm[0] * norm[0] + norm[1] * norm[1] + norm[2] * norm[2]);
509 // If point is on one edge, fn may be very small and the normal does not make sense - avoid divzero
510 if (fn < 1E-10) {
511 norm[0] = 1.;
512 norm[1] = 0.;
513 norm[2] = 0.;
514 } else {
515 norm[0] /= fn;
516 norm[1] /= fn;
517 norm[2] /= fn;
518 }
519 // Make sure the dot product of the normal and the direction is positive.
520 if (dir[0] > -2. && dir[0] * norm[0] + dir[1] * norm[1] + dir[2] * norm[2] < 0) {
521 norm[0] = -norm[0];
522 norm[1] = -norm[1];
523 norm[2] = -norm[2];
524 }
525}
526
527////////////////////////////////////////////////////////////////////////////////
528/// Test if point is inside this shape.
529
530Bool_t TGeoArb8::Contains(const Double_t *point) const
531{
532 // first check Z range
533 if (TMath::Abs(point[2]) > fDz)
534 return kFALSE;
535 // compute intersection between Z plane containing point and the arb8
536 Double_t poly[8];
537 // memset(&poly[0], 0, 8*sizeof(Double_t));
538 // SetPlaneVertices(point[2], &poly[0]);
539 Double_t cf = 0.5 * (fDz - point[2]) / fDz;
540 Int_t i;
541 for (i = 0; i < 4; i++) {
542 poly[2 * i] = fXY[i + 4][0] + cf * (fXY[i][0] - fXY[i + 4][0]);
543 poly[2 * i + 1] = fXY[i + 4][1] + cf * (fXY[i][1] - fXY[i + 4][1]);
544 }
545 return InsidePolygon(point[0], point[1], poly);
546}
547
548////////////////////////////////////////////////////////////////////////////////
549/// Computes distance to plane ipl :
550/// - ipl=0 : points 0,4,1,5
551/// - ipl=1 : points 1,5,2,6
552/// - ipl=2 : points 2,6,3,7
553/// - ipl=3 : points 3,7,0,4
554
555Double_t TGeoArb8::DistToPlane(const Double_t *point, const Double_t *dir, Int_t ipl, Bool_t in) const
556{
557 Double_t xa, xb, xc, xd;
558 Double_t ya, yb, yc, yd;
559 Double_t eps = 10. * TGeoShape::Tolerance();
560 Double_t norm[3];
561 Double_t dirp[3];
562 Double_t ndotd = 0;
563 Int_t j = (ipl + 1) % 4;
564 xa = fXY[ipl][0];
565 ya = fXY[ipl][1];
566 xb = fXY[ipl + 4][0];
567 yb = fXY[ipl + 4][1];
568 xc = fXY[j][0];
569 yc = fXY[j][1];
570 xd = fXY[4 + j][0];
571 yd = fXY[4 + j][1];
572 Double_t dz2 = 0.5 / fDz;
573 Double_t tx1 = dz2 * (xb - xa);
574 Double_t ty1 = dz2 * (yb - ya);
575 Double_t tx2 = dz2 * (xd - xc);
576 Double_t ty2 = dz2 * (yd - yc);
577 Double_t dzp = fDz + point[2];
578 Double_t xs1 = xa + tx1 * dzp;
579 Double_t ys1 = ya + ty1 * dzp;
580 Double_t xs2 = xc + tx2 * dzp;
581 Double_t ys2 = yc + ty2 * dzp;
582 Double_t dxs = xs2 - xs1;
583 Double_t dys = ys2 - ys1;
584 Double_t dtx = tx2 - tx1;
585 Double_t dty = ty2 - ty1;
586 Double_t a = (dtx * dir[1] - dty * dir[0] + (tx1 * ty2 - tx2 * ty1) * dir[2]) * dir[2];
587 Double_t signa = TMath::Sign(1., a);
588 Double_t b = dxs * dir[1] - dys * dir[0] +
589 (dtx * point[1] - dty * point[0] + ty2 * xs1 - ty1 * xs2 + tx1 * ys2 - tx2 * ys1) * dir[2];
590 Double_t c = dxs * point[1] - dys * point[0] + xs1 * ys2 - xs2 * ys1;
591 Double_t x1, x2, y1, y2, xp, yp, zi, s;
592 if (TMath::Abs(a) < eps) {
593 // Surface is planar
594 if (TMath::Abs(b) < eps)
595 return TGeoShape::Big(); // Track parallel to surface
596 s = -c / b;
597 if (TMath::Abs(s) < 1.E-6 && TMath::Abs(TMath::Abs(point[2]) - fDz) > eps) {
598 memcpy(dirp, dir, 3 * sizeof(Double_t));
599 dirp[0] = -3;
600 // Compute normal pointing outside
601 ((TGeoArb8 *)this)->ComputeNormal(point, dirp, norm);
602 ndotd = dir[0] * norm[0] + dir[1] * norm[1] + dir[2] * norm[2];
603 if (!in)
604 ndotd *= -1.;
605 if (ndotd > 0) {
606 s = TMath::Max(0., s);
607 zi = (point[0] - xs1) * (point[0] - xs2) + (point[1] - ys1) * (point[1] - ys2);
608 if (zi <= 0)
609 return s;
610 }
611 return TGeoShape::Big();
612 }
613 if (s < 0)
614 return TGeoShape::Big();
615 } else {
616 Double_t d = b * b - 4 * a * c;
617 if (d < 0)
618 return TGeoShape::Big();
619 Double_t smin = 0.5 * (-b - signa * TMath::Sqrt(d)) / a;
620 Double_t smax = 0.5 * (-b + signa * TMath::Sqrt(d)) / a;
621 s = smin;
622 if (TMath::Abs(s) < 1.E-6 && TMath::Abs(TMath::Abs(point[2]) - fDz) > eps) {
623 memcpy(dirp, dir, 3 * sizeof(Double_t));
624 dirp[0] = -3;
625 // Compute normal pointing outside
626 ((TGeoArb8 *)this)->ComputeNormal(point, dirp, norm);
627 ndotd = dir[0] * norm[0] + dir[1] * norm[1] + dir[2] * norm[2];
628 if (!in)
629 ndotd *= -1.;
630 if (ndotd > 0)
631 return TMath::Max(0., s);
632 s = 0.; // ignore
633 }
634 if (s > eps) {
635 // Check smin
636 zi = point[2] + s * dir[2];
637 if (TMath::Abs(zi) < fDz) {
638 x1 = xs1 + tx1 * dir[2] * s;
639 x2 = xs2 + tx2 * dir[2] * s;
640 xp = point[0] + s * dir[0];
641 y1 = ys1 + ty1 * dir[2] * s;
642 y2 = ys2 + ty2 * dir[2] * s;
643 yp = point[1] + s * dir[1];
644 zi = (xp - x1) * (xp - x2) + (yp - y1) * (yp - y2);
645 if (zi <= 0)
646 return s;
647 }
648 }
649 // Smin failed, try smax
650 s = smax;
651 if (TMath::Abs(s) < 1.E-6 && TMath::Abs(TMath::Abs(point[2]) - fDz) > eps) {
652 memcpy(dirp, dir, 3 * sizeof(Double_t));
653 dirp[0] = -3;
654 // Compute normal pointing outside
655 ((TGeoArb8 *)this)->ComputeNormal(point, dirp, norm);
656 ndotd = dir[0] * norm[0] + dir[1] * norm[1] + dir[2] * norm[2];
657 if (!in)
658 ndotd *= -1.;
659 if (ndotd > 0)
660 s = TMath::Max(0., s);
661 else
662 s = TGeoShape::Big();
663 return s;
664 }
665 }
666 if (s > eps) {
667 // Check smin
668 zi = point[2] + s * dir[2];
669 if (TMath::Abs(zi) < fDz) {
670 x1 = xs1 + tx1 * dir[2] * s;
671 x2 = xs2 + tx2 * dir[2] * s;
672 xp = point[0] + s * dir[0];
673 y1 = ys1 + ty1 * dir[2] * s;
674 y2 = ys2 + ty2 * dir[2] * s;
675 yp = point[1] + s * dir[1];
676 zi = (xp - x1) * (xp - x2) + (yp - y1) * (yp - y2);
677 if (zi <= 0)
678 return s;
679 }
680 }
681 return TGeoShape::Big();
682}
683
684////////////////////////////////////////////////////////////////////////////////
685/// Computes distance from outside point to surface of the shape.
686
687Double_t TGeoArb8::DistFromOutside(const Double_t *point, const Double_t *dir, Int_t /*iact*/, Double_t step,
688 Double_t * /*safe*/) const
689{
690 Double_t sdist = TGeoBBox::DistFromOutside(point, dir, fDX, fDY, fDZ, fOrigin, step);
691 if (sdist >= step)
692 return TGeoShape::Big();
693 Double_t snext;
694 // check Z planes
695 if (TMath::Abs(point[2]) > fDz - 1.E-8) {
696 Double_t pt[3];
697 if (point[2] * dir[2] < 0) {
698 pt[2] = fDz * TMath::Sign(1., point[2]);
699 snext = TMath::Max((pt[2] - point[2]) / dir[2], 0.);
700 for (Int_t j = 0; j < 2; j++)
701 pt[j] = point[j] + snext * dir[j];
702 if (Contains(&pt[0]))
703 return snext;
704 }
705 }
706 // check lateral faces
707 Double_t dist;
708 snext = TGeoShape::Big();
709 for (Int_t i = 0; i < 4; i++) {
710 dist = DistToPlane(point, dir, i, kFALSE);
711 if (dist < snext)
712 snext = dist;
713 }
714 return snext;
715}
716
717////////////////////////////////////////////////////////////////////////////////
718/// Compute distance from inside point to surface of the shape.
719
720Double_t TGeoArb8::DistFromInside(const Double_t *point, const Double_t *dir, Int_t /*iact*/, Double_t /*step*/,
721 Double_t * /*safe*/) const
722{
723 Int_t i;
724 Double_t distz = TGeoShape::Big();
725 Double_t distl = TGeoShape::Big();
726 Double_t dist;
727 Double_t pt[3] = {0., 0., 0.};
728 if (dir[2] < 0) {
729 distz = (-fDz - point[2]) / dir[2];
730 pt[2] = -fDz;
731 } else {
732 if (dir[2] > 0)
733 distz = (fDz - point[2]) / dir[2];
734 pt[2] = fDz;
735 }
736 for (i = 0; i < 4; i++) {
737 dist = DistToPlane(point, dir, i, kTRUE);
738 if (dist < distl)
739 distl = dist;
740 }
741 if (distz < distl) {
742 pt[0] = point[0] + distz * dir[0];
743 pt[1] = point[1] + distz * dir[1];
744 if (!Contains(pt))
745 distz = TGeoShape::Big();
746 }
747 dist = TMath::Min(distz, distl);
748 if (dist < 0 || dist > 1.E10)
749 return 0.;
750 return dist;
751#ifdef OLDALGORITHM
752 //#else
753 // compute distance to plane ipl :
754 // ipl=0 : points 0,4,1,5
755 // ipl=1 : points 1,5,2,6
756 // ipl=2 : points 2,6,3,7
757 // ipl=3 : points 3,7,0,4
758 Double_t distmin;
759 Bool_t lateral_cross = kFALSE;
760 if (dir[2] < 0) {
761 distmin = (-fDz - point[2]) / dir[2];
762 } else {
763 if (dir[2] > 0)
764 distmin = (fDz - point[2]) / dir[2];
765 else
766 distmin = TGeoShape::Big();
767 }
768 Double_t dz2 = 0.5 / fDz;
769 Double_t xa, xb, xc, xd;
770 Double_t ya, yb, yc, yd;
771 Double_t eps = 100. * TGeoShape::Tolerance();
772 for (Int_t ipl = 0; ipl < 4; ipl++) {
773 Int_t j = (ipl + 1) % 4;
774 xa = fXY[ipl][0];
775 ya = fXY[ipl][1];
776 xb = fXY[ipl + 4][0];
777 yb = fXY[ipl + 4][1];
778 xc = fXY[j][0];
779 yc = fXY[j][1];
780 xd = fXY[4 + j][0];
781 yd = fXY[4 + j][1];
782
783 Double_t tx1 = dz2 * (xb - xa);
784 Double_t ty1 = dz2 * (yb - ya);
785 Double_t tx2 = dz2 * (xd - xc);
786 Double_t ty2 = dz2 * (yd - yc);
787 Double_t dzp = fDz + point[2];
788 Double_t xs1 = xa + tx1 * dzp;
789 Double_t ys1 = ya + ty1 * dzp;
790 Double_t xs2 = xc + tx2 * dzp;
791 Double_t ys2 = yc + ty2 * dzp;
792 Double_t dxs = xs2 - xs1;
793 Double_t dys = ys2 - ys1;
794 Double_t dtx = tx2 - tx1;
795 Double_t dty = ty2 - ty1;
796 Double_t a = (dtx * dir[1] - dty * dir[0] + (tx1 * ty2 - tx2 * ty1) * dir[2]) * dir[2];
797 Double_t b = dxs * dir[1] - dys * dir[0] +
798 (dtx * point[1] - dty * point[0] + ty2 * xs1 - ty1 * xs2 + tx1 * ys2 - tx2 * ys1) * dir[2];
799 Double_t c = dxs * point[1] - dys * point[0] + xs1 * ys2 - xs2 * ys1;
800 Double_t s = TGeoShape::Big();
801 if (TMath::Abs(a) < eps) {
802 if (TMath::Abs(b) < eps)
803 continue;
804 s = -c / b;
805 if (s > eps && s < distmin) {
806 distmin = s;
807 lateral_cross = kTRUE;
808 }
809 continue;
810 }
811 Double_t d = b * b - 4 * a * c;
812 if (d >= 0.) {
813 if (a > 0)
814 s = 0.5 * (-b - TMath::Sqrt(d)) / a;
815 else
816 s = 0.5 * (-b + TMath::Sqrt(d)) / a;
817 if (s > eps) {
818 if (s < distmin) {
819 distmin = s;
820 lateral_cross = kTRUE;
821 }
822 } else {
823 if (a > 0)
824 s = 0.5 * (-b + TMath::Sqrt(d)) / a;
825 else
826 s = 0.5 * (-b - TMath::Sqrt(d)) / a;
827 if (s > eps && s < distmin) {
828 distmin = s;
829 lateral_cross = kTRUE;
830 }
831 }
832 }
833 }
834
835 if (!lateral_cross) {
836 // We have to make sure that track crosses the top or bottom.
837 if (distmin > 1.E10)
838 return TGeoShape::Tolerance();
839 Double_t pt[2];
840 pt[0] = point[0] + distmin * dir[0];
841 pt[1] = point[1] + distmin * dir[1];
842 // Check if propagated point is in the polygon
843 Double_t poly[8];
844 Int_t i = 0;
845 if (dir[2] > 0.)
846 i = 4;
847 for (Int_t j = 0; j < 4; j++) {
848 poly[2 * j] = fXY[j + i][0];
849 poly[2 * j + 1] = fXY[j + i][1];
850 }
851 if (!InsidePolygon(pt[0], pt[1], poly))
852 return TGeoShape::Tolerance();
853 }
854 return distmin;
855#endif
856}
857
858////////////////////////////////////////////////////////////////////////////////
859/// Divide this shape along one axis.
860
861TGeoVolume *TGeoArb8::Divide(TGeoVolume *voldiv, const char * /*divname*/, Int_t /*iaxis*/, Int_t /*ndiv*/,
862 Double_t /*start*/, Double_t /*step*/)
863{
864 Error("Divide", "Division of an arbitrary trapezoid not implemented");
865 return voldiv;
866}
867
868////////////////////////////////////////////////////////////////////////////////
869/// Get shape range on a given axis.
870
871Double_t TGeoArb8::GetAxisRange(Int_t iaxis, Double_t &xlo, Double_t &xhi) const
872{
873 xlo = 0;
874 xhi = 0;
875 Double_t dx = 0;
876 if (iaxis == 3) {
877 xlo = -fDz;
878 xhi = fDz;
879 dx = xhi - xlo;
880 return dx;
881 }
882 return dx;
883}
884
885////////////////////////////////////////////////////////////////////////////////
886/// Fill vector param[4] with the bounding cylinder parameters. The order
887/// is the following : Rmin, Rmax, Phi1, Phi2
888
889void TGeoArb8::GetBoundingCylinder(Double_t *param) const
890{
891 // first compute rmin/rmax
892 Double_t rmaxsq = 0;
893 Double_t rsq;
894 Int_t i;
895 for (i = 0; i < 8; i++) {
896 rsq = fXY[i][0] * fXY[i][0] + fXY[i][1] * fXY[i][1];
897 rmaxsq = TMath::Max(rsq, rmaxsq);
898 }
899 param[0] = 0.; // Rmin
900 param[1] = rmaxsq; // Rmax
901 param[2] = 0.; // Phi1
902 param[3] = 360.; // Phi2
903}
904
905////////////////////////////////////////////////////////////////////////////////
906/// Fills real parameters of a positioned box inside this arb8. Returns 0 if successful.
907
908Int_t TGeoArb8::GetFittingBox(const TGeoBBox *parambox, TGeoMatrix *mat, Double_t &dx, Double_t &dy, Double_t &dz) const
909{
910 dx = dy = dz = 0;
911 if (mat->IsRotation()) {
912 Error("GetFittingBox", "cannot handle parametrized rotated volumes");
913 return 1; // ### rotation not accepted ###
914 }
915 //--> translate the origin of the parametrized box to the frame of this box.
916 Double_t origin[3];
917 mat->LocalToMaster(parambox->GetOrigin(), origin);
918 if (!Contains(origin)) {
919 Error("GetFittingBox", "wrong matrix - parametrized box is outside this");
920 return 1; // ### wrong matrix ###
921 }
922 //--> now we have to get the valid range for all parametrized axis
923 Double_t dd[3];
924 dd[0] = parambox->GetDX();
925 dd[1] = parambox->GetDY();
926 dd[2] = parambox->GetDZ();
927 //-> check if Z range is fixed
928 if (dd[2] < 0) {
929 dd[2] = TMath::Min(origin[2] + fDz, fDz - origin[2]);
930 if (dd[2] < 0) {
931 Error("GetFittingBox", "wrong matrix");
932 return 1;
933 }
934 }
935 if (dd[0] >= 0 && dd[1] >= 0) {
936 dx = dd[0];
937 dy = dd[1];
938 dz = dd[2];
939 return 0;
940 }
941 //-> check now vertices at Z = origin[2] +/- dd[2]
942 Double_t upper[8];
943 Double_t lower[8];
944 SetPlaneVertices(origin[2] - dd[2], lower);
945 SetPlaneVertices(origin[2] + dd[2], upper);
946 for (Int_t iaxis = 0; iaxis < 2; iaxis++) {
947 if (dd[iaxis] >= 0)
948 continue;
949 Double_t ddmin = TGeoShape::Big();
950 for (Int_t ivert = 0; ivert < 4; ivert++) {
951 ddmin = TMath::Min(ddmin, TMath::Abs(origin[iaxis] - lower[2 * ivert + iaxis]));
952 ddmin = TMath::Min(ddmin, TMath::Abs(origin[iaxis] - upper[2 * ivert + iaxis]));
953 }
954 dd[iaxis] = ddmin;
955 }
956 dx = dd[0];
957 dy = dd[1];
958 dz = dd[2];
959 return 0;
960}
961
962////////////////////////////////////////////////////////////////////////////////
963/// Computes normal to plane defined by P1, P2 and P3
964
965void TGeoArb8::GetPlaneNormal(Double_t *p1, Double_t *p2, Double_t *p3, Double_t *norm)
966{
967 Double_t cross = 0.;
968 Double_t v1[3], v2[3];
969 Int_t i;
970 for (i = 0; i < 3; i++) {
971 v1[i] = p2[i] - p1[i];
972 v2[i] = p3[i] - p1[i];
973 }
974 norm[0] = v1[1] * v2[2] - v1[2] * v2[1];
975 cross += norm[0] * norm[0];
976 norm[1] = v1[2] * v2[0] - v1[0] * v2[2];
977 cross += norm[1] * norm[1];
978 norm[2] = v1[0] * v2[1] - v1[1] * v2[0];
979 cross += norm[2] * norm[2];
980 if (TMath::Abs(cross) < TGeoShape::Tolerance())
981 return;
982 cross = 1. / TMath::Sqrt(cross);
983 for (i = 0; i < 3; i++)
984 norm[i] *= cross;
985}
986
987////////////////////////////////////////////////////////////////////////////////
988/// Fills array with n random points located on the surface of indexed facet.
989/// The output array must be provided with a length of minimum 3*npoints. Returns
990/// true if operation succeeded.
991/// Possible index values:
992/// - 0 - all facets together
993/// - 1 to 6 - facet index from bottom to top Z
994
995Bool_t TGeoArb8::GetPointsOnFacet(Int_t /*index*/, Int_t /*npoints*/, Double_t * /* array */) const
996{
997 return kFALSE;
998 /*
999 if (index<0 || index>6) return kFALSE;
1000 if (index==0) {
1001 // Just generate same number of points on each facet
1002 Int_t npts = npoints/6.;
1003 Int_t count = 0;
1004 for (Int_t ifacet=0; ifacet<6; ifacet++) {
1005 if (GetPointsOnFacet(ifacet+1, npts, &array[3*count])) count += npts;
1006 if (ifacet<5) npts = (npoints-count)/(5.-ifacet);
1007 }
1008 if (count>0) return kTRUE;
1009 return kFALSE;
1010 }
1011 Double_t z, cf;
1012 Double_t xmin=TGeoShape::Big();
1013 Double_t xmax=-xmin;
1014 Double_t ymin=TGeoShape::Big();
1015 Double_t ymax=-ymin;
1016 Double_t dy=0.;
1017 Double_t poly[8];
1018 Double_t point[2];
1019 Int_t i;
1020 if (index==1 || index==6) {
1021 z = (index==1)?-fDz:fDz;
1022 cf = 0.5*(fDz-z)/fDz;
1023 for (i=0; i<4; i++) {
1024 poly[2*i] = fXY[i+4][0]+cf*(fXY[i][0]-fXY[i+4][0]);
1025 poly[2*i+1] = fXY[i+4][1]+cf*(fXY[i][1]-fXY[i+4][1]);
1026 xmin = TMath::Min(xmin, poly[2*i]);
1027 xmax = TMath::Max(xmax, poly[2*i]);
1028 ymin = TMath::Min(ymin, poly[2*i]);
1029 ymax = TMath::Max(ymax, poly[2*i]);
1030 }
1031 }
1032 Int_t nshoot = 0;
1033 Int_t nmiss = 0;
1034 for (i=0; i<npoints; i++) {
1035 Double_t *point = &array[3*i];
1036 switch (surfindex) {
1037 case 1:
1038 case 6:
1039 while (nmiss<1000) {
1040 point[0] = xmin + (xmax-xmin)*gRandom->Rndm();
1041 point[1] = ymin + (ymax-ymin)*gRandom->Rndm();
1042 }
1043
1044 return InsidePolygon(point[0],point[1],poly);
1045 */
1046}
1047
1048////////////////////////////////////////////////////////////////////////////////
1049/// Finds if a point in XY plane is inside the polygon defines by PTS.
1050
1051Bool_t TGeoArb8::InsidePolygon(Double_t x, Double_t y, Double_t *pts)
1052{
1053 Int_t i, j;
1054 Double_t x1, y1, x2, y2;
1055 Double_t cross;
1056 for (i = 0; i < 4; i++) {
1057 j = (i + 1) % 4;
1058 x1 = pts[i << 1];
1059 y1 = pts[(i << 1) + 1];
1060 x2 = pts[j << 1];
1061 y2 = pts[(j << 1) + 1];
1062 cross = (x - x1) * (y2 - y1) - (y - y1) * (x2 - x1);
1063 if (cross < 0)
1064 return kFALSE;
1065 }
1066 return kTRUE;
1067}
1068
1069////////////////////////////////////////////////////////////////////////////////
1070/// Prints shape parameters
1071
1072void TGeoArb8::InspectShape() const
1073{
1074 printf("*** Shape %s: TGeoArb8 ***\n", GetName());
1075 if (IsTwisted())
1076 printf(" = TWISTED\n");
1077 for (Int_t ip = 0; ip < 8; ip++) {
1078 printf(" point #%i : x=%11.5f y=%11.5f z=%11.5f\n", ip, fXY[ip][0], fXY[ip][1], fDz * ((ip < 4) ? -1 : 1));
1079 }
1080 printf(" Bounding box:\n");
1081 TGeoBBox::InspectShape();
1082}
1083
1084////////////////////////////////////////////////////////////////////////////////
1085/// Computes the closest distance from given point to this shape.
1086
1087Double_t TGeoArb8::Safety(const Double_t *point, Bool_t in) const
1088{
1089 Double_t safz = fDz - TMath::Abs(point[2]);
1090 if (!in)
1091 safz = -safz;
1092 Int_t iseg;
1093 Double_t safe = TGeoShape::Big();
1094 Double_t lsq, ssq, dx, dy, dpx, dpy, u;
1095 if (IsTwisted()) {
1096 if (!in) {
1097 if (!TGeoBBox::Contains(point))
1098 return TGeoBBox::Safety(point, kFALSE);
1099 }
1100 // Point is also in the bounding box ;-(
1101 // Compute closest distance to any segment
1102 Double_t vert[8];
1103 Double_t *p1, *p2;
1104 Int_t isegmin = 0;
1105 Double_t umin = 0.;
1106 SetPlaneVertices(point[2], vert);
1107 for (iseg = 0; iseg < 4; iseg++) {
1108 if (safe < TGeoShape::Tolerance())
1109 return 0.;
1110 p1 = &vert[2 * iseg];
1111 p2 = &vert[2 * ((iseg + 1) % 4)];
1112 dx = p2[0] - p1[0];
1113 dy = p2[1] - p1[1];
1114 dpx = point[0] - p1[0];
1115 dpy = point[1] - p1[1];
1116
1117 lsq = dx * dx + dy * dy;
1118 u = (dpx * dx + dpy * dy) / lsq;
1119 if (u > 1) {
1120 dpx = point[0] - p2[0];
1121 dpy = point[1] - p2[1];
1122 } else {
1123 if (u >= 0) {
1124 dpx -= u * dx;
1125 dpy -= u * dy;
1126 }
1127 }
1128 ssq = dpx * dpx + dpy * dpy;
1129 if (ssq < safe) {
1130 isegmin = iseg;
1131 umin = u;
1132 safe = ssq;
1133 }
1134 }
1135 if (umin < 0)
1136 umin = 0.;
1137 if (umin > 1) {
1138 isegmin = (isegmin + 1) % 4;
1139 umin = 0.;
1140 }
1141 Int_t i1 = isegmin;
1142 Int_t i2 = (isegmin + 1) % 4;
1143 Double_t dx1 = fXY[i2][0] - fXY[i1][0];
1144 Double_t dx2 = fXY[i2 + 4][0] - fXY[i1 + 4][0];
1145 Double_t dy1 = fXY[i2][1] - fXY[i1][1];
1146 Double_t dy2 = fXY[i2 + 4][1] - fXY[i1 + 4][1];
1147 dx = dx1 + umin * (dx2 - dx1);
1148 dy = dy1 + umin * (dy2 - dy1);
1149 safe *= 1. - 4. * fDz * fDz / (dx * dx + dy * dy + 4. * fDz * fDz);
1150 safe = TMath::Sqrt(safe);
1151 if (in)
1152 return TMath::Min(safz, safe);
1153 return TMath::Max(safz, safe);
1154 }
1155
1156 Double_t saf[5];
1157 saf[0] = safz;
1158
1159 for (iseg = 0; iseg < 4; iseg++)
1160 saf[iseg + 1] = SafetyToFace(point, iseg, in);
1161 if (in)
1162 safe = saf[TMath::LocMin(5, saf)];
1163 else
1164 safe = saf[TMath::LocMax(5, saf)];
1165 if (safe < 0)
1166 return 0.;
1167 return safe;
1168}
1169
1170////////////////////////////////////////////////////////////////////////////////
1171/// Estimate safety to lateral plane defined by segment iseg in range [0,3]
1172/// Might be negative: plane seen only from inside.
1173
1174Double_t TGeoArb8::SafetyToFace(const Double_t *point, Int_t iseg, Bool_t in) const
1175{
1176 Double_t vertices[12];
1177 Int_t ipln = (iseg + 1) % 4;
1178 // point 1
1179 vertices[0] = fXY[iseg][0];
1180 vertices[1] = fXY[iseg][1];
1181 vertices[2] = -fDz;
1182 // point 2
1183 vertices[3] = fXY[ipln][0];
1184 vertices[4] = fXY[ipln][1];
1185 vertices[5] = -fDz;
1186 // point 3
1187 vertices[6] = fXY[ipln + 4][0];
1188 vertices[7] = fXY[ipln + 4][1];
1189 vertices[8] = fDz;
1190 // point 4
1191 vertices[9] = fXY[iseg + 4][0];
1192 vertices[10] = fXY[iseg + 4][1];
1193 vertices[11] = fDz;
1194 Double_t safe;
1195 Double_t norm[3];
1196 Double_t *p1, *p2, *p3;
1197 p1 = &vertices[0];
1198 p2 = &vertices[9];
1199 p3 = &vertices[6];
1200 if (IsSamePoint(p2, p3)) {
1201 p3 = &vertices[3];
1202 if (IsSamePoint(p1, p3))
1203 return -TGeoShape::Big(); // skip single segment
1204 }
1205 GetPlaneNormal(p1, p2, p3, norm);
1206 safe = (point[0] - p1[0]) * norm[0] + (point[1] - p1[1]) * norm[1] + (point[2] - p1[2]) * norm[2];
1207 if (in)
1208 return (-safe);
1209 return safe;
1210}
1211
1212////////////////////////////////////////////////////////////////////////////////
1213/// Save a primitive as a C++ statement(s) on output stream "out".
1214
1215void TGeoArb8::SavePrimitive(std::ostream &out, Option_t * /*option*/ /*= ""*/)
1216{
1217 if (TObject::TestBit(kGeoSavePrimitive))
1218 return;
1219 out << " // Shape: " << GetName() << " type: " << ClassName() << std::endl;
1220 out << " dz = " << fDz << ";" << std::endl;
1221 out << " vert[0] = " << fXY[0][0] << ";" << std::endl;
1222 out << " vert[1] = " << fXY[0][1] << ";" << std::endl;
1223 out << " vert[2] = " << fXY[1][0] << ";" << std::endl;
1224 out << " vert[3] = " << fXY[1][1] << ";" << std::endl;
1225 out << " vert[4] = " << fXY[2][0] << ";" << std::endl;
1226 out << " vert[5] = " << fXY[2][1] << ";" << std::endl;
1227 out << " vert[6] = " << fXY[3][0] << ";" << std::endl;
1228 out << " vert[7] = " << fXY[3][1] << ";" << std::endl;
1229 out << " vert[8] = " << fXY[4][0] << ";" << std::endl;
1230 out << " vert[9] = " << fXY[4][1] << ";" << std::endl;
1231 out << " vert[10] = " << fXY[5][0] << ";" << std::endl;
1232 out << " vert[11] = " << fXY[5][1] << ";" << std::endl;
1233 out << " vert[12] = " << fXY[6][0] << ";" << std::endl;
1234 out << " vert[13] = " << fXY[6][1] << ";" << std::endl;
1235 out << " vert[14] = " << fXY[7][0] << ";" << std::endl;
1236 out << " vert[15] = " << fXY[7][1] << ";" << std::endl;
1237 out << " TGeoShape *" << GetPointerName() << " = new TGeoArb8(\"" << GetName() << "\", dz,vert);" << std::endl;
1238 TObject::SetBit(TGeoShape::kGeoSavePrimitive);
1239}
1240
1241////////////////////////////////////////////////////////////////////////////////
1242/// Computes intersection points between plane at zpl and non-horizontal edges.
1243
1244void TGeoArb8::SetPlaneVertices(Double_t zpl, Double_t *vertices) const
1245{
1246 Double_t cf = 0.5 * (fDz - zpl) / fDz;
1247 for (Int_t i = 0; i < 4; i++) {
1248 vertices[2 * i] = fXY[i + 4][0] + cf * (fXY[i][0] - fXY[i + 4][0]);
1249 vertices[2 * i + 1] = fXY[i + 4][1] + cf * (fXY[i][1] - fXY[i + 4][1]);
1250 }
1251}
1252
1253////////////////////////////////////////////////////////////////////////////////
1254/// Set all arb8 params in one step.
1255/// param[0] = dz
1256/// param[1] = x0
1257/// param[2] = y0
1258/// ...
1259
1260void TGeoArb8::SetDimensions(Double_t *param)
1261{
1262 fDz = param[0];
1263 for (Int_t i = 0; i < 8; i++) {
1264 fXY[i][0] = param[2 * i + 1];
1265 fXY[i][1] = param[2 * i + 2];
1266 }
1267 ComputeTwist();
1268 ComputeBBox();
1269}
1270
1271////////////////////////////////////////////////////////////////////////////////
1272/// Creates arb8 mesh points
1273
1274void TGeoArb8::SetPoints(Double_t *points) const
1275{
1276 for (Int_t i = 0; i < 8; i++) {
1277 points[3 * i] = fXY[i][0];
1278 points[3 * i + 1] = fXY[i][1];
1279 points[3 * i + 2] = (i < 4) ? -fDz : fDz;
1280 }
1281}
1282
1283////////////////////////////////////////////////////////////////////////////////
1284/// Creates arb8 mesh points
1285
1286void TGeoArb8::SetPoints(Float_t *points) const
1287{
1288 for (Int_t i = 0; i < 8; i++) {
1289 points[3 * i] = fXY[i][0];
1290 points[3 * i + 1] = fXY[i][1];
1291 points[3 * i + 2] = (i < 4) ? -fDz : fDz;
1292 }
1293}
1294
1295////////////////////////////////////////////////////////////////////////////////
1296/// Set values for a given vertex.
1297
1298void TGeoArb8::SetVertex(Int_t vnum, Double_t x, Double_t y)
1299{
1300 if (vnum < 0 || vnum > 7) {
1301 Error("SetVertex", "Invalid vertex number");
1302 return;
1303 }
1304 fXY[vnum][0] = x;
1305 fXY[vnum][1] = y;
1306 if (vnum == 7) {
1307 ComputeTwist();
1308 ComputeBBox();
1309 }
1310}
1311
1312////////////////////////////////////////////////////////////////////////////////
1313/// Fill size of this 3-D object
1314
1315void TGeoArb8::Sizeof3D() const
1316{
1317 TGeoBBox::Sizeof3D();
1318}
1319
1320////////////////////////////////////////////////////////////////////////////////
1321/// Stream an object of class TGeoManager.
1322
1323void TGeoArb8::Streamer(TBuffer &R__b)
1324{
1325 if (R__b.IsReading()) {
1326 R__b.ReadClassBuffer(TGeoArb8::Class(), this);
1327 ComputeTwist();
1328 } else {
1329 R__b.WriteClassBuffer(TGeoArb8::Class(), this);
1330 }
1331}
1332
1333////////////////////////////////////////////////////////////////////////////////
1334/// Check the inside status for each of the points in the array.
1335/// Input: Array of point coordinates + vector size
1336/// Output: Array of Booleans for the inside of each point
1337
1338void TGeoArb8::Contains_v(const Double_t *points, Bool_t *inside, Int_t vecsize) const
1339{
1340 for (Int_t i = 0; i < vecsize; i++)
1341 inside[i] = Contains(&points[3 * i]);
1342}
1343
1344////////////////////////////////////////////////////////////////////////////////
1345/// Compute the normal for an array o points so that norm.dot.dir is positive
1346/// Input: Arrays of point coordinates and directions + vector size
1347/// Output: Array of normal directions
1348
1349void TGeoArb8::ComputeNormal_v(const Double_t *points, const Double_t *dirs, Double_t *norms, Int_t vecsize)
1350{
1351 for (Int_t i = 0; i < vecsize; i++)
1352 ComputeNormal(&points[3 * i], &dirs[3 * i], &norms[3 * i]);
1353}
1354
1355////////////////////////////////////////////////////////////////////////////////
1356/// Compute distance from array of input points having directions specified by dirs. Store output in dists
1357
1358void TGeoArb8::DistFromInside_v(const Double_t *points, const Double_t *dirs, Double_t *dists, Int_t vecsize,
1359 Double_t *step) const
1360{
1361 for (Int_t i = 0; i < vecsize; i++)
1362 dists[i] = DistFromInside(&points[3 * i], &dirs[3 * i], 3, step[i]);
1363}
1364
1365////////////////////////////////////////////////////////////////////////////////
1366/// Compute distance from array of input points having directions specified by dirs. Store output in dists
1367
1368void TGeoArb8::DistFromOutside_v(const Double_t *points, const Double_t *dirs, Double_t *dists, Int_t vecsize,
1369 Double_t *step) const
1370{
1371 for (Int_t i = 0; i < vecsize; i++)
1372 dists[i] = DistFromOutside(&points[3 * i], &dirs[3 * i], 3, step[i]);
1373}
1374
1375////////////////////////////////////////////////////////////////////////////////
1376/// Compute safe distance from each of the points in the input array.
1377/// Input: Array of point coordinates, array of statuses for these points, size of the arrays
1378/// Output: Safety values
1379
1380void TGeoArb8::Safety_v(const Double_t *points, const Bool_t *inside, Double_t *safe, Int_t vecsize) const
1381{
1382 for (Int_t i = 0; i < vecsize; i++)
1383 safe[i] = Safety(&points[3 * i], inside[i]);
1384}
1385
1386
1387////////////////////////////////////////////////////////////////////////////////
1388/// Default ctor
1389
1390TGeoTrap::TGeoTrap()
1391{
1392 fDz = 0;
1393 fTheta = 0;
1394 fPhi = 0;
1395 fH1 = fH2 = fBl1 = fBl2 = fTl1 = fTl2 = fAlpha1 = fAlpha2 = 0;
1396}
1397
1398////////////////////////////////////////////////////////////////////////////////
1399/// Constructor providing just a range in Z, theta and phi.
1400
1401TGeoTrap::TGeoTrap(Double_t dz, Double_t theta, Double_t phi) : TGeoArb8("", 0, nullptr)
1402{
1403 fDz = dz;
1404 fTheta = theta;
1405 fPhi = phi;
1406 fH1 = fH2 = fBl1 = fBl2 = fTl1 = fTl2 = fAlpha1 = fAlpha2 = 0;
1407}
1408
1409////////////////////////////////////////////////////////////////////////////////
1410/// Normal constructor.
1411
1412TGeoTrap::TGeoTrap(Double_t dz, Double_t theta, Double_t phi, Double_t h1, Double_t bl1, Double_t tl1, Double_t alpha1,
1413 Double_t h2, Double_t bl2, Double_t tl2, Double_t alpha2)
1414 : TGeoArb8("", 0, nullptr)
1415{
1416 fDz = dz;
1417 fTheta = theta;
1418 fPhi = phi;
1419 fH1 = h1;
1420 fH2 = h2;
1421 fBl1 = bl1;
1422 fBl2 = bl2;
1423 fTl1 = tl1;
1424 fTl2 = tl2;
1425 fAlpha1 = alpha1;
1426 fAlpha2 = alpha2;
1427 Double_t tx = TMath::Tan(theta * TMath::DegToRad()) * TMath::Cos(phi * TMath::DegToRad());
1428 Double_t ty = TMath::Tan(theta * TMath::DegToRad()) * TMath::Sin(phi * TMath::DegToRad());
1429 Double_t ta1 = TMath::Tan(alpha1 * TMath::DegToRad());
1430 Double_t ta2 = TMath::Tan(alpha2 * TMath::DegToRad());
1431 fXY[0][0] = -dz * tx - h1 * ta1 - bl1;
1432 fXY[0][1] = -dz * ty - h1;
1433 fXY[1][0] = -dz * tx + h1 * ta1 - tl1;
1434 fXY[1][1] = -dz * ty + h1;
1435 fXY[2][0] = -dz * tx + h1 * ta1 + tl1;
1436 fXY[2][1] = -dz * ty + h1;
1437 fXY[3][0] = -dz * tx - h1 * ta1 + bl1;
1438 fXY[3][1] = -dz * ty - h1;
1439 fXY[4][0] = dz * tx - h2 * ta2 - bl2;
1440 fXY[4][1] = dz * ty - h2;
1441 fXY[5][0] = dz * tx + h2 * ta2 - tl2;
1442 fXY[5][1] = dz * ty + h2;
1443 fXY[6][0] = dz * tx + h2 * ta2 + tl2;
1444 fXY[6][1] = dz * ty + h2;
1445 fXY[7][0] = dz * tx - h2 * ta2 + bl2;
1446 fXY[7][1] = dz * ty - h2;
1447 ComputeTwist();
1448 if ((dz < 0) || (h1 < 0) || (bl1 < 0) || (tl1 < 0) || (h2 < 0) || (bl2 < 0) || (tl2 < 0)) {
1449 SetShapeBit(kGeoRunTimeShape);
1450 } else
1451 TGeoArb8::ComputeBBox();
1452}
1453
1454////////////////////////////////////////////////////////////////////////////////
1455/// Constructor with name.
1456
1457TGeoTrap::TGeoTrap(const char *name, Double_t dz, Double_t theta, Double_t phi, Double_t h1, Double_t bl1, Double_t tl1,
1458 Double_t alpha1, Double_t h2, Double_t bl2, Double_t tl2, Double_t alpha2)
1459 : TGeoArb8(name, 0, nullptr)
1460{
1461 SetName(name);
1462 fDz = dz;
1463 fTheta = theta;
1464 fPhi = phi;
1465 fH1 = h1;
1466 fH2 = h2;
1467 fBl1 = bl1;
1468 fBl2 = bl2;
1469 fTl1 = tl1;
1470 fTl2 = tl2;
1471 fAlpha1 = alpha1;
1472 fAlpha2 = alpha2;
1473 for (Int_t i = 0; i < 8; i++) {
1474 fXY[i][0] = 0.0;
1475 fXY[i][1] = 0.0;
1476 }
1477 Double_t tx = TMath::Tan(theta * TMath::DegToRad()) * TMath::Cos(phi * TMath::DegToRad());
1478 Double_t ty = TMath::Tan(theta * TMath::DegToRad()) * TMath::Sin(phi * TMath::DegToRad());
1479 Double_t ta1 = TMath::Tan(alpha1 * TMath::DegToRad());
1480 Double_t ta2 = TMath::Tan(alpha2 * TMath::DegToRad());
1481 fXY[0][0] = -dz * tx - h1 * ta1 - bl1;
1482 fXY[0][1] = -dz * ty - h1;
1483 fXY[1][0] = -dz * tx + h1 * ta1 - tl1;
1484 fXY[1][1] = -dz * ty + h1;
1485 fXY[2][0] = -dz * tx + h1 * ta1 + tl1;
1486 fXY[2][1] = -dz * ty + h1;
1487 fXY[3][0] = -dz * tx - h1 * ta1 + bl1;
1488 fXY[3][1] = -dz * ty - h1;
1489 fXY[4][0] = dz * tx - h2 * ta2 - bl2;
1490 fXY[4][1] = dz * ty - h2;
1491 fXY[5][0] = dz * tx + h2 * ta2 - tl2;
1492 fXY[5][1] = dz * ty + h2;
1493 fXY[6][0] = dz * tx + h2 * ta2 + tl2;
1494 fXY[6][1] = dz * ty + h2;
1495 fXY[7][0] = dz * tx - h2 * ta2 + bl2;
1496 fXY[7][1] = dz * ty - h2;
1497 ComputeTwist();
1498 if ((dz < 0) || (h1 < 0) || (bl1 < 0) || (tl1 < 0) || (h2 < 0) || (bl2 < 0) || (tl2 < 0)) {
1499 SetShapeBit(kGeoRunTimeShape);
1500 } else
1501 TGeoArb8::ComputeBBox();
1502}
1503
1504////////////////////////////////////////////////////////////////////////////////
1505/// Destructor.
1506
1507TGeoTrap::~TGeoTrap() {}
1508
1509////////////////////////////////////////////////////////////////////////////////
1510/// Compute distance from inside point to surface of the trapezoid
1511
1512Double_t
1513TGeoTrap::DistFromInside(const Double_t *point, const Double_t *dir, Int_t iact, Double_t step, Double_t *safe) const
1514{
1515 if (iact < 3 && safe) {
1516 // compute safe distance
1517 *safe = Safety(point, kTRUE);
1518 if (iact == 0)
1519 return TGeoShape::Big();
1520 if (iact == 1 && step < *safe)
1521 return TGeoShape::Big();
1522 }
1523 // compute distance to get outside this shape
1524 // return TGeoArb8::DistFromInside(point, dir, iact, step, safe);
1525
1526 // compute distance to plane ipl :
1527 // ipl=0 : points 0,4,1,5
1528 // ipl=1 : points 1,5,2,6
1529 // ipl=2 : points 2,6,3,7
1530 // ipl=3 : points 3,7,0,4
1531 Double_t distmin;
1532 if (dir[2] < 0) {
1533 distmin = (-fDz - point[2]) / dir[2];
1534 } else {
1535 if (dir[2] > 0)
1536 distmin = (fDz - point[2]) / dir[2];
1537 else
1538 distmin = TGeoShape::Big();
1539 }
1540 Double_t xa, xb, xc;
1541 Double_t ya, yb, yc;
1542 for (Int_t ipl = 0; ipl < 4; ipl++) {
1543 Int_t j = (ipl + 1) % 4;
1544 xa = fXY[ipl][0];
1545 ya = fXY[ipl][1];
1546 xb = fXY[ipl + 4][0];
1547 yb = fXY[ipl + 4][1];
1548 xc = fXY[j][0];
1549 yc = fXY[j][1];
1550 Double_t ax, ay, az;
1551 ax = xb - xa;
1552 ay = yb - ya;
1553 az = 2. * fDz;
1554 Double_t bx, by;
1555 bx = xc - xa;
1556 by = yc - ya;
1557 Double_t ddotn = -dir[0] * az * by + dir[1] * az * bx + dir[2] * (ax * by - ay * bx);
1558 if (ddotn <= 0)
1559 continue; // entering
1560 Double_t saf = -(point[0] - xa) * az * by + (point[1] - ya) * az * bx + (point[2] + fDz) * (ax * by - ay * bx);
1561 if (saf >= 0.0)
1562 return 0.0;
1563 Double_t s = -saf / ddotn;
1564 if (s < distmin)
1565 distmin = s;
1566 }
1567 return distmin;
1568}
1569
1570////////////////////////////////////////////////////////////////////////////////
1571/// Compute distance from outside point to surface of the trapezoid
1572
1573Double_t
1574TGeoTrap::DistFromOutside(const Double_t *point, const Double_t *dir, Int_t iact, Double_t step, Double_t *safe) const
1575{
1576 if (iact < 3 && safe) {
1577 // compute safe distance
1578 *safe = Safety(point, kFALSE);
1579 if (iact == 0)
1580 return TGeoShape::Big();
1581 if (iact == 1 && step < *safe)
1582 return TGeoShape::Big();
1583 }
1584 // Check if the bounding box is crossed within the requested distance
1585 Double_t sdist = TGeoBBox::DistFromOutside(point, dir, fDX, fDY, fDZ, fOrigin, step);
1586 if (sdist >= step)
1587 return TGeoShape::Big();
1588 // compute distance to get outside this shape
1589 Bool_t in = kTRUE;
1590 Double_t pts[8];
1591 Double_t xnew, ynew, znew;
1592 Int_t i, j;
1593 if (point[2] < -fDz + TGeoShape::Tolerance()) {
1594 if (dir[2] < TGeoShape::Tolerance())
1595 return TGeoShape::Big();
1596 in = kFALSE;
1597 Double_t snxt = -(fDz + point[2]) / dir[2];
1598 xnew = point[0] + snxt * dir[0];
1599 ynew = point[1] + snxt * dir[1];
1600 for (i = 0; i < 4; i++) {
1601 j = i << 1;
1602 pts[j] = fXY[i][0];
1603 pts[j + 1] = fXY[i][1];
1604 }
1605 if (InsidePolygon(xnew, ynew, pts))
1606 return snxt;
1607 } else if (point[2] > fDz - TGeoShape::Tolerance()) {
1608 if (dir[2] > -TGeoShape::Tolerance())
1609 return TGeoShape::Big();
1610 in = kFALSE;
1611 Double_t snxt = (fDz - point[2]) / dir[2];
1612 xnew = point[0] + snxt * dir[0];
1613 ynew = point[1] + snxt * dir[1];
1614 for (i = 0; i < 4; i++) {
1615 j = i << 1;
1616 pts[j] = fXY[i + 4][0];
1617 pts[j + 1] = fXY[i + 4][1];
1618 }
1619 if (InsidePolygon(xnew, ynew, pts))
1620 return snxt;
1621 }
1622 // check lateral faces
1623 Double_t dz2 = 0.5 / fDz;
1624 Double_t xa, xb, xc, xd;
1625 Double_t ya, yb, yc, yd;
1626 Double_t ax, ay, az;
1627 Double_t bx, by;
1628 Double_t ddotn, saf;
1629 Double_t safmin = TGeoShape::Big();
1630 Bool_t exiting = kFALSE;
1631 for (i = 0; i < 4; i++) {
1632 j = (i + 1) % 4;
1633 xa = fXY[i][0];
1634 ya = fXY[i][1];
1635 xb = fXY[i + 4][0];
1636 yb = fXY[i + 4][1];
1637 xc = fXY[j][0];
1638 yc = fXY[j][1];
1639 xd = fXY[4 + j][0];
1640 yd = fXY[4 + j][1];
1641 ax = xb - xa;
1642 ay = yb - ya;
1643 az = 2. * fDz;
1644 bx = xc - xa;
1645 by = yc - ya;
1646 ddotn = -dir[0] * az * by + dir[1] * az * bx + dir[2] * (ax * by - ay * bx);
1647 saf = (point[0] - xa) * az * by - (point[1] - ya) * az * bx - (point[2] + fDz) * (ax * by - ay * bx);
1648
1649 if (saf <= 0) {
1650 // face visible from point outside
1651 in = kFALSE;
1652 if (ddotn >= 0)
1653 return TGeoShape::Big();
1654 Double_t snxt = saf / ddotn;
1655 znew = point[2] + snxt * dir[2];
1656 if (TMath::Abs(znew) <= fDz) {
1657 xnew = point[0] + snxt * dir[0];
1658 ynew = point[1] + snxt * dir[1];
1659 Double_t tx1 = dz2 * (xb - xa);
1660 Double_t ty1 = dz2 * (yb - ya);
1661 Double_t tx2 = dz2 * (xd - xc);
1662 Double_t ty2 = dz2 * (yd - yc);
1663 Double_t dzp = fDz + znew;
1664 Double_t xs1 = xa + tx1 * dzp;
1665 Double_t ys1 = ya + ty1 * dzp;
1666 Double_t xs2 = xc + tx2 * dzp;
1667 Double_t ys2 = yc + ty2 * dzp;
1668 if (TMath::Abs(xs1 - xs2) > TMath::Abs(ys1 - ys2)) {
1669 if ((xnew - xs1) * (xs2 - xnew) >= 0)
1670 return snxt;
1671 } else {
1672 if ((ynew - ys1) * (ys2 - ynew) >= 0)
1673 return snxt;
1674 }
1675 }
1676 } else {
1677 if (saf < safmin) {
1678 safmin = saf;
1679 if (ddotn >= 0)
1680 exiting = kTRUE;
1681 else
1682 exiting = kFALSE;
1683 }
1684 }
1685 }
1686 // Check also Z boundaries (point may be inside and close to Z) - Christian Hammann
1687 saf = fDz - TMath::Abs(point[2]);
1688 if (saf > 0 && saf < safmin)
1689 exiting = (point[2] * dir[2] > 0) ? kTRUE : kFALSE;
1690 if (!in)
1691 return TGeoShape::Big();
1692 if (exiting)
1693 return TGeoShape::Big();
1694 return 0.0;
1695}
1696
1697////////////////////////////////////////////////////////////////////////////////
1698/// Divide this trapezoid shape belonging to volume "voldiv" into ndiv volumes
1699/// called divname, from start position with the given step. Only Z divisions
1700/// are supported. For Z divisions just return the pointer to the volume to be
1701/// divided. In case a wrong division axis is supplied, returns pointer to
1702/// volume that was divided.
1703
1704TGeoVolume *
1705TGeoTrap::Divide(TGeoVolume *voldiv, const char *divname, Int_t iaxis, Int_t ndiv, Double_t start, Double_t step)
1706{
1707 TGeoShape *shape; //--- shape to be created
1708 TGeoVolume *vol; //--- division volume to be created
1709 TGeoVolumeMulti *vmulti; //--- generic divided volume
1710 TGeoPatternFinder *finder; //--- finder to be attached
1711 TString opt = ""; //--- option to be attached
1712 if (iaxis != 3) {
1713 Error("Divide", "cannot divide trapezoids on other axis than Z");
1714 return nullptr;
1715 }
1716 Double_t end = start + ndiv * step;
1717 Double_t points_lo[8];
1718 Double_t points_hi[8];
1719 finder = new TGeoPatternTrapZ(voldiv, ndiv, start, end);
1720 voldiv->SetFinder(finder);
1721 finder->SetDivIndex(voldiv->GetNdaughters());
1722 opt = "Z";
1723 vmulti = gGeoManager->MakeVolumeMulti(divname, voldiv->GetMedium());
1724 Double_t txz = ((TGeoPatternTrapZ *)finder)->GetTxz();
1725 Double_t tyz = ((TGeoPatternTrapZ *)finder)->GetTyz();
1726 Double_t zmin, zmax, ox, oy, oz;
1727 for (Int_t idiv = 0; idiv < ndiv; idiv++) {
1728 zmin = start + idiv * step;
1729 zmax = start + (idiv + 1) * step;
1730 oz = start + idiv * step + step / 2;
1731 ox = oz * txz;
1732 oy = oz * tyz;
1733 SetPlaneVertices(zmin, &points_lo[0]);
1734 SetPlaneVertices(zmax, &points_hi[0]);
1735 shape = new TGeoTrap(step / 2, fTheta, fPhi);
1736 for (Int_t vert1 = 0; vert1 < 4; vert1++)
1737 ((TGeoArb8 *)shape)->SetVertex(vert1, points_lo[2 * vert1] - ox, points_lo[2 * vert1 + 1] - oy);
1738 for (Int_t vert2 = 0; vert2 < 4; vert2++)
1739 ((TGeoArb8 *)shape)->SetVertex(vert2 + 4, points_hi[2 * vert2] - ox, points_hi[2 * vert2 + 1] - oy);
1740 vol = new TGeoVolume(divname, shape, voldiv->GetMedium());
1741 vmulti->AddVolume(vol);
1742 voldiv->AddNodeOffset(vol, idiv, oz, opt.Data());
1743 ((TGeoNodeOffset *)voldiv->GetNodes()->At(voldiv->GetNdaughters() - 1))->SetFinder(finder);
1744 }
1745 return vmulti;
1746}
1747
1748////////////////////////////////////////////////////////////////////////////////
1749/// In case shape has some negative parameters, these have to be computed
1750/// in order to fit the mother.
1751
1752TGeoShape *TGeoTrap::GetMakeRuntimeShape(TGeoShape *mother, TGeoMatrix * /*mat*/) const
1753{
1754 if (!TestShapeBit(kGeoRunTimeShape))
1755 return nullptr;
1756 if (mother->IsRunTimeShape()) {
1757 Error("GetMakeRuntimeShape", "invalid mother");
1758 return nullptr;
1759 }
1760 Double_t dz, h1, bl1, tl1, h2, bl2, tl2;
1761 if (fDz < 0)
1762 dz = ((TGeoTrap *)mother)->GetDz();
1763 else
1764 dz = fDz;
1765
1766 if (fH1 < 0)
1767 h1 = ((TGeoTrap *)mother)->GetH1();
1768 else
1769 h1 = fH1;
1770
1771 if (fH2 < 0)
1772 h2 = ((TGeoTrap *)mother)->GetH2();
1773 else
1774 h2 = fH2;
1775
1776 if (fBl1 < 0)
1777 bl1 = ((TGeoTrap *)mother)->GetBl1();
1778 else
1779 bl1 = fBl1;
1780
1781 if (fBl2 < 0)
1782 bl2 = ((TGeoTrap *)mother)->GetBl2();
1783 else
1784 bl2 = fBl2;
1785
1786 if (fTl1 < 0)
1787 tl1 = ((TGeoTrap *)mother)->GetTl1();
1788 else
1789 tl1 = fTl1;
1790
1791 if (fTl2 < 0)
1792 tl2 = ((TGeoTrap *)mother)->GetTl2();
1793 else
1794 tl2 = fTl2;
1795
1796 return (new TGeoTrap(dz, fTheta, fPhi, h1, bl1, tl1, fAlpha1, h2, bl2, tl2, fAlpha2));
1797}
1798
1799////////////////////////////////////////////////////////////////////////////////
1800/// Computes the closest distance from given point to this shape.
1801
1802Double_t TGeoTrap::Safety(const Double_t *point, Bool_t in) const
1803{
1804 Double_t saf[5];
1805 Double_t norm[3]; // normal to current facette
1806 Int_t i, j; // current facette index
1807 Double_t x0, y0, z0 = -fDz, x1, y1, z1 = fDz, x2, y2;
1808 Double_t ax, ay, az = z1 - z0, bx, by;
1809 Double_t fn, safe;
1810 //---> compute safety for lateral planes
1811 for (i = 0; i < 4; i++) {
1812 if (in)
1813 saf[i] = TGeoShape::Big();
1814 else
1815 saf[i] = 0.;
1816 x0 = fXY[i][0];
1817 y0 = fXY[i][1];
1818 x1 = fXY[i + 4][0];
1819 y1 = fXY[i + 4][1];
1820 ax = x1 - x0;
1821 ay = y1 - y0;
1822 j = (i + 1) % 4;
1823 x2 = fXY[j][0];
1824 y2 = fXY[j][1];
1825 bx = x2 - x0;
1826 by = y2 - y0;
1827 if (TMath::Abs(bx) < TGeoShape::Tolerance() && TMath::Abs(by) < TGeoShape::Tolerance()) {
1828 x2 = fXY[4 + j][0];
1829 y2 = fXY[4 + j][1];
1830 bx = x2 - x1;
1831 by = y2 - y1;
1832 if (TMath::Abs(bx) < TGeoShape::Tolerance() && TMath::Abs(by) < TGeoShape::Tolerance())
1833 continue;
1834 }
1835 norm[0] = -az * by;
1836 norm[1] = az * bx;
1837 norm[2] = ax * by - ay * bx;
1838 fn = TMath::Sqrt(norm[0] * norm[0] + norm[1] * norm[1] + norm[2] * norm[2]);
1839 if (fn < 1E-10)
1840 continue;
1841 saf[i] = (x0 - point[0]) * norm[0] + (y0 - point[1]) * norm[1] + (-fDz - point[2]) * norm[2];
1842 if (in) {
1843 saf[i] = TMath::Abs(saf[i]) / fn; // they should be all positive anyway
1844 } else {
1845 saf[i] = -saf[i] / fn; // only negative values are interesting
1846 }
1847 }
1848 saf[4] = fDz - TMath::Abs(point[2]);
1849 if (in) {
1850 safe = saf[0];
1851 for (j = 1; j < 5; j++)
1852 if (saf[j] < safe)
1853 safe = saf[j];
1854 } else {
1855 saf[4] = -saf[4];
1856 safe = saf[0];
1857 for (j = 1; j < 5; j++)
1858 if (saf[j] > safe)
1859 safe = saf[j];
1860 }
1861 return safe;
1862}
1863
1864////////////////////////////////////////////////////////////////////////////////
1865/// Save a primitive as a C++ statement(s) on output stream "out".
1866
1867void TGeoTrap::SavePrimitive(std::ostream &out, Option_t * /*option*/ /*= ""*/)
1868{
1869 if (TObject::TestBit(kGeoSavePrimitive))
1870 return;
1871 out << " // Shape: " << GetName() << " type: " << ClassName() << std::endl;
1872 out << " dz = " << fDz << ";" << std::endl;
1873 out << " theta = " << fTheta << ";" << std::endl;
1874 out << " phi = " << fPhi << ";" << std::endl;
1875 out << " h1 = " << fH1 << ";" << std::endl;
1876 out << " bl1 = " << fBl1 << ";" << std::endl;
1877 out << " tl1 = " << fTl1 << ";" << std::endl;
1878 out << " alpha1 = " << fAlpha1 << ";" << std::endl;
1879 out << " h2 = " << fH2 << ";" << std::endl;
1880 out << " bl2 = " << fBl2 << ";" << std::endl;
1881 out << " tl2 = " << fTl2 << ";" << std::endl;
1882 out << " alpha2 = " << fAlpha2 << ";" << std::endl;
1883 out << " TGeoShape *" << GetPointerName() << " = new TGeoTrap(\"" << GetName()
1884 << "\", dz,theta,phi,h1,bl1,tl1,alpha1,h2,bl2,tl2,alpha2);" << std::endl;
1885 TObject::SetBit(TGeoShape::kGeoSavePrimitive);
1886}
1887
1888////////////////////////////////////////////////////////////////////////////////
1889/// Set all arb8 params in one step.
1890/// - param[0] = dz
1891/// - param[1] = theta
1892/// - param[2] = phi
1893/// - param[3] = h1
1894/// - param[4] = bl1
1895/// - param[5] = tl1
1896/// - param[6] = alpha1
1897/// - param[7] = h2
1898/// - param[8] = bl2
1899/// - param[9] = tl2
1900/// - param[10] = alpha2
1901
1902void TGeoTrap::SetDimensions(Double_t *param)
1903{
1904 fDz = param[0];
1905 fTheta = param[1];
1906 fPhi = param[2];
1907 fH1 = param[3];
1908 fH2 = param[7];
1909 fBl1 = param[4];
1910 fBl2 = param[8];
1911 fTl1 = param[5];
1912 fTl2 = param[9];
1913 fAlpha1 = param[6];
1914 fAlpha2 = param[10];
1915 Double_t tx = TMath::Tan(fTheta * TMath::DegToRad()) * TMath::Cos(fPhi * TMath::DegToRad());
1916 Double_t ty = TMath::Tan(fTheta * TMath::DegToRad()) * TMath::Sin(fPhi * TMath::DegToRad());
1917 Double_t ta1 = TMath::Tan(fAlpha1 * TMath::DegToRad());
1918 Double_t ta2 = TMath::Tan(fAlpha2 * TMath::DegToRad());
1919 fXY[0][0] = -fDz * tx - fH1 * ta1 - fBl1;
1920 fXY[0][1] = -fDz * ty - fH1;
1921 fXY[1][0] = -fDz * tx + fH1 * ta1 - fTl1;
1922 fXY[1][1] = -fDz * ty + fH1;
1923 fXY[2][0] = -fDz * tx + fH1 * ta1 + fTl1;
1924 fXY[2][1] = -fDz * ty + fH1;
1925 fXY[3][0] = -fDz * tx - fH1 * ta1 + fBl1;
1926 fXY[3][1] = -fDz * ty - fH1;
1927 fXY[4][0] = fDz * tx - fH2 * ta2 - fBl2;
1928 fXY[4][1] = fDz * ty - fH2;
1929 fXY[5][0] = fDz * tx + fH2 * ta2 - fTl2;
1930 fXY[5][1] = fDz * ty + fH2;
1931 fXY[6][0] = fDz * tx + fH2 * ta2 + fTl2;
1932 fXY[6][1] = fDz * ty + fH2;
1933 fXY[7][0] = fDz * tx - fH2 * ta2 + fBl2;
1934 fXY[7][1] = fDz * ty - fH2;
1935 ComputeTwist();
1936 if ((fDz < 0) || (fH1 < 0) || (fBl1 < 0) || (fTl1 < 0) || (fH2 < 0) || (fBl2 < 0) || (fTl2 < 0)) {
1937 SetShapeBit(kGeoRunTimeShape);
1938 } else
1939 TGeoArb8::ComputeBBox();
1940}
1941
1942////////////////////////////////////////////////////////////////////////////////
1943/// Compute distance from array of input points having directions specified by dirs. Store output in dists
1944
1945void TGeoTrap::DistFromInside_v(const Double_t *points, const Double_t *dirs, Double_t *dists, Int_t vecsize,
1946 Double_t *step) const
1947{
1948 for (Int_t i = 0; i < vecsize; i++)
1949 dists[i] = DistFromInside(&points[3 * i], &dirs[3 * i], 3, step[i]);
1950}
1951
1952////////////////////////////////////////////////////////////////////////////////
1953/// Compute distance from array of input points having directions specified by dirs. Store output in dists
1954
1955void TGeoTrap::DistFromOutside_v(const Double_t *points, const Double_t *dirs, Double_t *dists, Int_t vecsize,
1956 Double_t *step) const
1957{
1958 for (Int_t i = 0; i < vecsize; i++)
1959 dists[i] = DistFromOutside(&points[3 * i], &dirs[3 * i], 3, step[i]);
1960}
1961
1962////////////////////////////////////////////////////////////////////////////////
1963/// Compute safe distance from each of the points in the input array.
1964/// Input: Array of point coordinates, array of statuses for these points, size of the arrays
1965/// Output: Safety values
1966
1967void TGeoTrap::Safety_v(const Double_t *points, const Bool_t *inside, Double_t *safe, Int_t vecsize) const
1968{
1969 for (Int_t i = 0; i < vecsize; i++)
1970 safe[i] = Safety(&points[3 * i], inside[i]);
1971}
1972
1973
1974////////////////////////////////////////////////////////////////////////////////
1975/// Default ctor
1976
1977TGeoGtra::TGeoGtra()
1978{
1979 fTwistAngle = 0;
1980}
1981
1982////////////////////////////////////////////////////////////////////////////////
1983/// Constructor.
1984
1985TGeoGtra::TGeoGtra(Double_t dz, Double_t theta, Double_t phi, Double_t twist, Double_t h1, Double_t bl1, Double_t tl1,
1986 Double_t alpha1, Double_t h2, Double_t bl2, Double_t tl2, Double_t alpha2)
1987 : TGeoTrap(dz, theta, phi, h1, bl1, tl1, alpha1, h2, bl2, tl2, alpha2)
1988{
1989 fTwistAngle = twist;
1990 Double_t x, y;
1991 Double_t th = theta * TMath::DegToRad();
1992 Double_t ph = phi * TMath::DegToRad();
1993 // Coordinates of the center of the bottom face
1994 Double_t xc = -dz * TMath::Sin(th) * TMath::Cos(ph);
1995 Double_t yc = -dz * TMath::Sin(th) * TMath::Sin(ph);
1996
1997 Int_t i;
1998
1999 for (i = 0; i < 4; i++) {
2000 x = fXY[i][0] - xc;
2001 y = fXY[i][1] - yc;
2002 fXY[i][0] =
2003 x * TMath::Cos(-0.5 * twist * TMath::DegToRad()) + y * TMath::Sin(-0.5 * twist * TMath::DegToRad()) + xc;
2004 fXY[i][1] =
2005 -x * TMath::Sin(-0.5 * twist * TMath::DegToRad()) + y * TMath::Cos(-0.5 * twist * TMath::DegToRad()) + yc;
2006 }
2007 // Coordinates of the center of the top face
2008 xc = -xc;
2009 yc = -yc;
2010 for (i = 4; i < 8; i++) {
2011 x = fXY[i][0] - xc;
2012 y = fXY[i][1] - yc;
2013 fXY[i][0] =
2014 x * TMath::Cos(0.5 * twist * TMath::DegToRad()) + y * TMath::Sin(0.5 * twist * TMath::DegToRad()) + xc;
2015 fXY[i][1] =
2016 -x * TMath::Sin(0.5 * twist * TMath::DegToRad()) + y * TMath::Cos(0.5 * twist * TMath::DegToRad()) + yc;
2017 }
2018 ComputeTwist();
2019 if ((dz < 0) || (h1 < 0) || (bl1 < 0) || (tl1 < 0) || (h2 < 0) || (bl2 < 0) || (tl2 < 0))
2020 SetShapeBit(kGeoRunTimeShape);
2021 else
2022 TGeoArb8::ComputeBBox();
2023}
2024
2025////////////////////////////////////////////////////////////////////////////////
2026/// Constructor providing the name of the shape.
2027
2028TGeoGtra::TGeoGtra(const char *name, Double_t dz, Double_t theta, Double_t phi, Double_t twist, Double_t h1,
2029 Double_t bl1, Double_t tl1, Double_t alpha1, Double_t h2, Double_t bl2, Double_t tl2,
2030 Double_t alpha2)
2031 : TGeoTrap(name, dz, theta, phi, h1, bl1, tl1, alpha1, h2, bl2, tl2, alpha2)
2032{
2033 fTwistAngle = twist;
2034 Double_t x, y;
2035 Double_t th = theta * TMath::DegToRad();
2036 Double_t ph = phi * TMath::DegToRad();
2037 // Coordinates of the center of the bottom face
2038 Double_t xc = -dz * TMath::Sin(th) * TMath::Cos(ph);
2039 Double_t yc = -dz * TMath::Sin(th) * TMath::Sin(ph);
2040
2041 Int_t i;
2042
2043 for (i = 0; i < 4; i++) {
2044 x = fXY[i][0] - xc;
2045 y = fXY[i][1] - yc;
2046 fXY[i][0] =
2047 x * TMath::Cos(-0.5 * twist * TMath::DegToRad()) + y * TMath::Sin(-0.5 * twist * TMath::DegToRad()) + xc;
2048 fXY[i][1] =
2049 -x * TMath::Sin(-0.5 * twist * TMath::DegToRad()) + y * TMath::Cos(-0.5 * twist * TMath::DegToRad()) + yc;
2050 }
2051 // Coordinates of the center of the top face
2052 xc = -xc;
2053 yc = -yc;
2054 for (i = 4; i < 8; i++) {
2055 x = fXY[i][0] - xc;
2056 y = fXY[i][1] - yc;
2057 fXY[i][0] =
2058 x * TMath::Cos(0.5 * twist * TMath::DegToRad()) + y * TMath::Sin(0.5 * twist * TMath::DegToRad()) + xc;
2059 fXY[i][1] =
2060 -x * TMath::Sin(0.5 * twist * TMath::DegToRad()) + y * TMath::Cos(0.5 * twist * TMath::DegToRad()) + yc;
2061 }
2062 ComputeTwist();
2063 if ((dz < 0) || (h1 < 0) || (bl1 < 0) || (tl1 < 0) || (h2 < 0) || (bl2 < 0) || (tl2 < 0))
2064 SetShapeBit(kGeoRunTimeShape);
2065 else
2066 TGeoArb8::ComputeBBox();
2067}
2068
2069////////////////////////////////////////////////////////////////////////////////
2070/// Destructor.
2071
2072TGeoGtra::~TGeoGtra() {}
2073
2074////////////////////////////////////////////////////////////////////////////////
2075/// Compute distance from inside point to surface of the shape.
2076
2077Double_t
2078TGeoGtra::DistFromInside(const Double_t *point, const Double_t *dir, Int_t iact, Double_t step, Double_t *safe) const
2079{
2080 if (iact < 3 && safe) {
2081 // compute safe distance
2082 *safe = Safety(point, kTRUE);
2083 if (iact == 0)
2084 return TGeoShape::Big();
2085 if (iact == 1 && step < *safe)
2086 return TGeoShape::Big();
2087 }
2088 // compute distance to get outside this shape
2089 return TGeoArb8::DistFromInside(point, dir, iact, step, safe);
2090}
2091
2092////////////////////////////////////////////////////////////////////////////////
2093/// Compute distance from inside point to surface of the shape.
2094
2095Double_t
2096TGeoGtra::DistFromOutside(const Double_t *point, const Double_t *dir, Int_t iact, Double_t step, Double_t *safe) const
2097{
2098 if (iact < 3 && safe) {
2099 // compute safe distance
2100 *safe = Safety(point, kTRUE);
2101 if (iact == 0)
2102 return TGeoShape::Big();
2103 if (iact == 1 && step < *safe)
2104 return TGeoShape::Big();
2105 }
2106 // compute distance to get outside this shape
2107 return TGeoArb8::DistFromOutside(point, dir, iact, step, safe);
2108}
2109
2110////////////////////////////////////////////////////////////////////////////////
2111/// In case shape has some negative parameters, these has to be computed
2112/// in order to fit the mother
2113
2114TGeoShape *TGeoGtra::GetMakeRuntimeShape(TGeoShape *mother, TGeoMatrix * /*mat*/) const
2115{
2116 if (!TestShapeBit(kGeoRunTimeShape))
2117 return nullptr;
2118 if (mother->IsRunTimeShape()) {
2119 Error("GetMakeRuntimeShape", "invalid mother");
2120 return nullptr;
2121 }
2122 Double_t dz, h1, bl1, tl1, h2, bl2, tl2;
2123 if (fDz < 0)
2124 dz = ((TGeoTrap *)mother)->GetDz();
2125 else
2126 dz = fDz;
2127 if (fH1 < 0)
2128 h1 = ((TGeoTrap *)mother)->GetH1();
2129 else
2130 h1 = fH1;
2131 if (fH2 < 0)
2132 h2 = ((TGeoTrap *)mother)->GetH2();
2133 else
2134 h2 = fH2;
2135 if (fBl1 < 0)
2136 bl1 = ((TGeoTrap *)mother)->GetBl1();
2137 else
2138 bl1 = fBl1;
2139 if (fBl2 < 0)
2140 bl2 = ((TGeoTrap *)mother)->GetBl2();
2141 else
2142 bl2 = fBl2;
2143 if (fTl1 < 0)
2144 tl1 = ((TGeoTrap *)mother)->GetTl1();
2145 else
2146 tl1 = fTl1;
2147 if (fTl2 < 0)
2148 tl2 = ((TGeoTrap *)mother)->GetTl2();
2149 else
2150 tl2 = fTl2;
2151 return (new TGeoGtra(dz, fTheta, fPhi, fTwistAngle, h1, bl1, tl1, fAlpha1, h2, bl2, tl2, fAlpha2));
2152}
2153
2154////////////////////////////////////////////////////////////////////////////////
2155/// Computes the closest distance from given point to this shape.
2156
2157Double_t TGeoGtra::Safety(const Double_t *point, Bool_t in) const
2158{
2159 return TGeoArb8::Safety(point, in);
2160}
2161
2162////////////////////////////////////////////////////////////////////////////////
2163/// Save a primitive as a C++ statement(s) on output stream "out".
2164
2165void TGeoGtra::SavePrimitive(std::ostream &out, Option_t * /*option*/ /*= ""*/)
2166{
2167 if (TObject::TestBit(kGeoSavePrimitive))
2168 return;
2169 out << " // Shape: " << GetName() << " type: " << ClassName() << std::endl;
2170 out << " dz = " << fDz << ";" << std::endl;
2171 out << " theta = " << fTheta << ";" << std::endl;
2172 out << " phi = " << fPhi << ";" << std::endl;
2173 out << " twist = " << fTwistAngle << ";" << std::endl;
2174 out << " h1 = " << fH1 << ";" << std::endl;
2175 out << " bl1 = " << fBl1 << ";" << std::endl;
2176 out << " tl1 = " << fTl1 << ";" << std::endl;
2177 out << " alpha1 = " << fAlpha1 << ";" << std::endl;
2178 out << " h2 = " << fH2 << ";" << std::endl;
2179 out << " bl2 = " << fBl2 << ";" << std::endl;
2180 out << " tl2 = " << fTl2 << ";" << std::endl;
2181 out << " alpha2 = " << fAlpha2 << ";" << std::endl;
2182 out << " TGeoShape *" << GetPointerName() << " = new TGeoGtra(\"" << GetName()
2183 << "\", dz,theta,phi,twist,h1,bl1,tl1,alpha1,h2,bl2,tl2,alpha2);" << std::endl;
2184 TObject::SetBit(TGeoShape::kGeoSavePrimitive);
2185}
2186
2187////////////////////////////////////////////////////////////////////////////////
2188/// Set all arb8 params in one step.
2189/// - param[0] = dz
2190/// - param[1] = theta
2191/// - param[2] = phi
2192/// - param[3] = h1
2193/// - param[4] = bl1
2194/// - param[5] = tl1
2195/// - param[6] = alpha1
2196/// - param[7] = h2
2197/// - param[8] = bl2
2198/// - param[9] = tl2
2199/// - param[10] = alpha2
2200/// - param[11] = twist
2201
2202void TGeoGtra::SetDimensions(Double_t *param)
2203{
2204 TGeoTrap::SetDimensions(param);
2205 fTwistAngle = param[11];
2206 Double_t x, y;
2207 Double_t twist = fTwistAngle;
2208 Double_t th = fTheta * TMath::DegToRad();
2209 Double_t ph = fPhi * TMath::DegToRad();
2210 // Coordinates of the center of the bottom face
2211 Double_t xc = -fDz * TMath::Sin(th) * TMath::Cos(ph);
2212 Double_t yc = -fDz * TMath::Sin(th) * TMath::Sin(ph);
2213
2214 Int_t i;
2215
2216 for (i = 0; i < 4; i++) {
2217 x = fXY[i][0] - xc;
2218 y = fXY[i][1] - yc;
2219 fXY[i][0] =
2220 x * TMath::Cos(-0.5 * twist * TMath::DegToRad()) + y * TMath::Sin(-0.5 * twist * TMath::DegToRad()) + xc;
2221 fXY[i][1] =
2222 -x * TMath::Sin(-0.5 * twist * TMath::DegToRad()) + y * TMath::Cos(-0.5 * twist * TMath::DegToRad()) + yc;
2223 }
2224 // Coordinates of the center of the top face
2225 xc = -xc;
2226 yc = -yc;
2227 for (i = 4; i < 8; i++) {
2228 x = fXY[i][0] - xc;
2229 y = fXY[i][1] - yc;
2230 fXY[i][0] =
2231 x * TMath::Cos(0.5 * twist * TMath::DegToRad()) + y * TMath::Sin(0.5 * twist * TMath::DegToRad()) + xc;
2232 fXY[i][1] =
2233 -x * TMath::Sin(0.5 * twist * TMath::DegToRad()) + y * TMath::Cos(0.5 * twist * TMath::DegToRad()) + yc;
2234 }
2235 ComputeTwist();
2236 if ((fDz < 0) || (fH1 < 0) || (fBl1 < 0) || (fTl1 < 0) || (fH2 < 0) || (fBl2 < 0) || (fTl2 < 0))
2237 SetShapeBit(kGeoRunTimeShape);
2238 else
2239 TGeoArb8::ComputeBBox();
2240}
2241
2242////////////////////////////////////////////////////////////////////////////////
2243/// Compute distance from array of input points having directions specified by dirs. Store output in dists
2244
2245void TGeoGtra::DistFromInside_v(const Double_t *points, const Double_t *dirs, Double_t *dists, Int_t vecsize,
2246 Double_t *step) const
2247{
2248 for (Int_t i = 0; i < vecsize; i++)
2249 dists[i] = DistFromInside(&points[3 * i], &dirs[3 * i], 3, step[i]);
2250}
2251
2252////////////////////////////////////////////////////////////////////////////////
2253/// Compute distance from array of input points having directions specified by dirs. Store output in dists
2254
2255void TGeoGtra::DistFromOutside_v(const Double_t *points, const Double_t *dirs, Double_t *dists, Int_t vecsize,
2256 Double_t *step) const
2257{
2258 for (Int_t i = 0; i < vecsize; i++)
2259 dists[i] = DistFromOutside(&points[3 * i], &dirs[3 * i], 3, step[i]);
2260}
2261
2262////////////////////////////////////////////////////////////////////////////////
2263/// Compute safe distance from each of the points in the input array.
2264/// Input: Array of point coordinates, array of statuses for these points, size of the arrays
2265/// Output: Safety values
2266
2267void TGeoGtra::Safety_v(const Double_t *points, const Bool_t *inside, Double_t *safe, Int_t vecsize) const
2268{
2269 for (Int_t i = 0; i < vecsize; i++)
2270 safe[i] = Safety(&points[3 * i], inside[i]);
2271}
d
#define d(i)
Definition RSha256.hxx:102
b
#define b(i)
Definition RSha256.hxx:100
c
#define c(i)
Definition RSha256.hxx:101
a
#define a(i)
Definition RSha256.hxx:99
start
start
Definition Rotated.cxx:223
Int_t
int Int_t
Signed integer 4 bytes (int).
Definition RtypesCore.h:59
Bool_t
bool Bool_t
Boolean (0=false, 1=true) (bool).
Definition RtypesCore.h:77
kFALSE
constexpr Bool_t kFALSE
Definition RtypesCore.h:108
Double_t
double Double_t
Double 8 bytes.
Definition RtypesCore.h:73
Float_t
float Float_t
Float 4 bytes (float).
Definition RtypesCore.h:71
kTRUE
constexpr Bool_t kTRUE
Definition RtypesCore.h:107
Option_t
const char Option_t
Option string (const char).
Definition RtypesCore.h:80
Error
Error("WriteTObject","The current directory (%s) is not associated with a file. The object (%s) has not been written.", GetName(), objname)
name
char name[80]
Definition TGX11.cxx:148
gGeoManager
externTGeoManager * gGeoManager
Definition TGeoManager.h:620
xmin
float xmin
Definition THbookFile.cxx:95
ymin
float ymin
Definition THbookFile.cxx:95
xmax
float xmax
Definition THbookFile.cxx:95
ymax
float ymax
Definition THbookFile.cxx:95
points
point * points
Definition X3DBuffer.c:22
TBuffer
Buffer base class used for serializing objects.
Definition TBuffer.h:43
TBuffer::ReadClassBuffer
virtual Int_t ReadClassBuffer(const TClass *cl, void *pointer, const TClass *onfile_class=nullptr)=0
TBuffer::IsReading
Bool_t IsReading() const
Definition TBuffer.h:86
TBuffer::WriteClassBuffer
virtual Int_t WriteClassBuffer(const TClass *cl, void *pointer)=0
TGeoArb8::Safety_v
void Safety_v(const Double_t *points, const Bool_t *inside, Double_t *safe, Int_t vecsize) const override
TGeoArb8::ComputeNormal
void ComputeNormal(const Double_t *point, const Double_t *dir, Double_t *norm) const override
TGeoArb8::DistFromOutside_v
void DistFromOutside_v(const Double_t *points, const Double_t *dirs, Double_t *dists, Int_t vecsize, Double_t *step) const override
TGeoArb8::GetClosestEdge
Double_t GetClosestEdge(const Double_t *point, Double_t *vert, Int_t &isegment) const
TGeoArb8::fXY
Double_t fXY[8][2]
Definition TGeoArb8.h:30
TGeoArb8::SavePrimitive
void SavePrimitive(std::ostream &out, Option_t *option="") override
Save a primitive as a C++ statement(s) on output stream "out".
TGeoArb8::SetPoints
void SetPoints(Double_t *points) const override
TGeoArb8::SetVertex
virtual void SetVertex(Int_t vnum, Double_t x, Double_t y)
TGeoArb8::CopyTwist
void CopyTwist(Double_t *twist=nullptr)
TGeoArb8::SetDimensions
void SetDimensions(Double_t *param) override
TGeoArb8::Contains
Bool_t Contains(const Double_t *point) const override
TGeoArb8::GetPointsOnFacet
Bool_t GetPointsOnFacet(Int_t, Int_t, Double_t *) const override
TGeoArb8::TGeoArb8
TGeoArb8(const TGeoArb8 &)=delete
TGeoArb8::ComputeTwist
void ComputeTwist()
TGeoArb8::DistToPlane
Double_t DistToPlane(const Double_t *point, const Double_t *dir, Int_t ipl, Bool_t in) const
TGeoArb8::Streamer
void Streamer(TBuffer &) override
Stream an object of class TObject.
TGeoArb8::InsidePolygon
static Bool_t InsidePolygon(Double_t x, Double_t y, Double_t *pts)
TGeoArb8::Divide
TGeoVolume * Divide(TGeoVolume *voldiv, const char *divname, Int_t iaxis, Int_t ndiv, Double_t start, Double_t step) override
TGeoArb8::fDz
Double_t fDz
Definition TGeoArb8.h:28
TGeoArb8::GetBoundingCylinder
void GetBoundingCylinder(Double_t *param) const override
TGeoArb8::fTwist
Double_t * fTwist
! [4] tangents of twist angles
Definition TGeoArb8.h:29
TGeoArb8::GetAxisRange
Double_t GetAxisRange(Int_t iaxis, Double_t &xlo, Double_t &xhi) const override
TGeoArb8::GetPlaneNormal
static void GetPlaneNormal(Double_t *p1, Double_t *p2, Double_t *p3, Double_t *norm)
TGeoArb8::Capacity
Double_t Capacity() const override
TGeoArb8::SafetyToFace
Double_t SafetyToFace(const Double_t *point, Int_t iseg, Bool_t in) const
TGeoArb8::Sizeof3D
void Sizeof3D() const override
TGeoArb8::DistFromOutside
Double_t DistFromOutside(const Double_t *point, const Double_t *dir, Int_t iact=1, Double_t step=TGeoShape::Big(), Double_t *safe=nullptr) const override
TGeoArb8::IsTwisted
Bool_t IsTwisted() const
Definition TGeoArb8.h:83
TGeoArb8::Contains_v
void Contains_v(const Double_t *points, Bool_t *inside, Int_t vecsize) const override
TGeoArb8::InspectShape
void InspectShape() const override
TGeoArb8::Safety
Double_t Safety(const Double_t *point, Bool_t in=kTRUE) const override
TGeoArb8::GetTwist
Double_t GetTwist(Int_t iseg) const
TGeoArb8::DistFromInside
Double_t DistFromInside(const Double_t *point, const Double_t *dir, Int_t iact=1, Double_t step=TGeoShape::Big(), Double_t *safe=nullptr) const override
TGeoArb8::ComputeNormal_v
void ComputeNormal_v(const Double_t *points, const Double_t *dirs, Double_t *norms, Int_t vecsize) override
TGeoArb8::IsSamePoint
static Bool_t IsSamePoint(const Double_t *p1, const Double_t *p2)
Definition TGeoArb8.h:77
TGeoArb8::~TGeoArb8
~TGeoArb8() override
TGeoArb8::DistFromInside_v
void DistFromInside_v(const Double_t *points, const Double_t *dirs, Double_t *dists, Int_t vecsize, Double_t *step) const override
TGeoArb8::Class
static TClass * Class()
TGeoArb8::SetPlaneVertices
void SetPlaneVertices(Double_t zpl, Double_t *vertices) const
TGeoArb8::ComputeBBox
void ComputeBBox() override
TGeoArb8::GetFittingBox
Int_t GetFittingBox(const TGeoBBox *parambox, TGeoMatrix *mat, Double_t &dx, Double_t &dy, Double_t &dz) const override
TGeoBBox::GetOrigin
virtual const Double_t * GetOrigin() const
Definition TGeoBBox.h:81
TGeoBBox::fDX
Double_t fDX
Definition TGeoBBox.h:21
TGeoBBox::DistFromOutside
Double_t DistFromOutside(const Double_t *point, const Double_t *dir, Int_t iact=1, Double_t step=TGeoShape::Big(), Double_t *safe=nullptr) const override
TGeoBBox::GetDX
virtual Double_t GetDX() const
Definition TGeoBBox.h:78
TGeoBBox::GetDZ
virtual Double_t GetDZ() const
Definition TGeoBBox.h:80
TGeoBBox::GetDY
virtual Double_t GetDY() const
Definition TGeoBBox.h:79
TGeoBBox::fOrigin
Double_t fOrigin[3]
Definition TGeoBBox.h:24
TGeoBBox::InspectShape
void InspectShape() const override
TGeoBBox::Contains
Bool_t Contains(const Double_t *point) const override
TGeoBBox::Safety
Double_t Safety(const Double_t *point, Bool_t in=kTRUE) const override
TGeoBBox::Sizeof3D
void Sizeof3D() const override
TGeoBBox::fDY
Double_t fDY
Definition TGeoBBox.h:22
TGeoBBox::fDZ
Double_t fDZ
Definition TGeoBBox.h:23
TGeoGtra::SetDimensions
void SetDimensions(Double_t *param) override
TGeoGtra::DistFromInside
Double_t DistFromInside(const Double_t *point, const Double_t *dir, Int_t iact=1, Double_t step=TGeoShape::Big(), Double_t *safe=nullptr) const override
TGeoGtra::fTwistAngle
Double_t fTwistAngle
Definition TGeoArb8.h:156
TGeoGtra::SavePrimitive
void SavePrimitive(std::ostream &out, Option_t *option="") override
Save a primitive as a C++ statement(s) on output stream "out".
TGeoGtra::DistFromInside_v
void DistFromInside_v(const Double_t *points, const Double_t *dirs, Double_t *dists, Int_t vecsize, Double_t *step) const override
TGeoGtra::GetMakeRuntimeShape
TGeoShape * GetMakeRuntimeShape(TGeoShape *mother, TGeoMatrix *mat) const override
TGeoGtra::Safety
Double_t Safety(const Double_t *point, Bool_t in=kTRUE) const override
TGeoGtra::Safety_v
void Safety_v(const Double_t *points, const Bool_t *inside, Double_t *safe, Int_t vecsize) const override
TGeoGtra::DistFromOutside_v
void DistFromOutside_v(const Double_t *points, const Double_t *dirs, Double_t *dists, Int_t vecsize, Double_t *step) const override
TGeoGtra::DistFromOutside
Double_t DistFromOutside(const Double_t *point, const Double_t *dir, Int_t iact=1, Double_t step=TGeoShape::Big(), Double_t *safe=nullptr) const override
TGeoGtra::~TGeoGtra
~TGeoGtra() override
TGeoMatrix
Geometrical transformation package.
Definition TGeoMatrix.h:39
TGeoMatrix::IsRotation
Bool_t IsRotation() const
Definition TGeoMatrix.h:66
TGeoMatrix::LocalToMaster
virtual void LocalToMaster(const Double_t *local, Double_t *master) const
convert a point by multiplying its column vector (x, y, z, 1) to matrix inverse
Definition TGeoMatrix.cxx:366
TGeoPatternFinder::SetDivIndex
void SetDivIndex(Int_t index)
Definition TGeoPatternFinder.h:88
TGeoShape
Base abstract class for all shapes.
Definition TGeoShape.h:25
TGeoShape::Big
static Double_t Big()
Definition TGeoShape.h:95
TGeoShape::IsSegCrossing
static Bool_t IsSegCrossing(Double_t x1, Double_t y1, Double_t x2, Double_t y2, Double_t x3, Double_t y3, Double_t x4, Double_t y4)
Check if segments (A,B) and (C,D) are crossing, where: A(x1,y1), B(x2,y2), C(x3,y3),...
Definition TGeoShape.cxx:355
TGeoShape::SetShapeBit
void SetShapeBit(UInt_t f, Bool_t set)
Equivalent of TObject::SetBit.
Definition TGeoShape.cxx:568
TGeoShape::TGeoShape
TGeoShape()
Default constructor.
Definition TGeoShape.cxx:165
TGeoShape::IsSameWithinTolerance
static Bool_t IsSameWithinTolerance(Double_t a, Double_t b)
Check if two numbers differ with less than a tolerance.
Definition TGeoShape.cxx:344
TGeoShape::IsRunTimeShape
Bool_t IsRunTimeShape() const
Definition TGeoShape.h:152
TGeoShape::GetPointerName
const char * GetPointerName() const
Provide a pointer name containing uid.
Definition TGeoShape.cxx:752
TGeoShape::GetName
const char * GetName() const override
Get the shape name.
Definition TGeoShape.cxx:248
TGeoShape::kGeoClosedShape
@ kGeoClosedShape
Definition TGeoShape.h:59
TGeoShape::kGeoSavePrimitive
@ kGeoSavePrimitive
Definition TGeoShape.h:65
TGeoShape::kGeoRunTimeShape
@ kGeoRunTimeShape
Definition TGeoShape.h:40
TGeoShape::Tolerance
static Double_t Tolerance()
Definition TGeoShape.h:98
TGeoShape::TestShapeBit
Bool_t TestShapeBit(UInt_t f) const
Definition TGeoShape.h:177
TGeoTrap::~TGeoTrap
~TGeoTrap() override
TGeoTrap::fTl2
Double_t fTl2
Definition TGeoArb8.h:110
TGeoTrap::fBl2
Double_t fBl2
Definition TGeoArb8.h:109
TGeoTrap::fAlpha2
Double_t fAlpha2
Definition TGeoArb8.h:111
TGeoTrap::SetDimensions
void SetDimensions(Double_t *param) override
TGeoTrap::DistFromInside
Double_t DistFromInside(const Double_t *point, const Double_t *dir, Int_t iact=1, Double_t step=TGeoShape::Big(), Double_t *safe=nullptr) const override
TGeoTrap::Safety
Double_t Safety(const Double_t *point, Bool_t in=kTRUE) const override
TGeoTrap::SavePrimitive
void SavePrimitive(std::ostream &out, Option_t *option="") override
Save a primitive as a C++ statement(s) on output stream "out".
TGeoTrap::fBl1
Double_t fBl1
Definition TGeoArb8.h:105
TGeoTrap::DistFromInside_v
void DistFromInside_v(const Double_t *points, const Double_t *dirs, Double_t *dists, Int_t vecsize, Double_t *step) const override
TGeoTrap::fH1
Double_t fH1
Definition TGeoArb8.h:104
TGeoTrap::Safety_v
void Safety_v(const Double_t *points, const Bool_t *inside, Double_t *safe, Int_t vecsize) const override
TGeoTrap::fPhi
Double_t fPhi
Definition TGeoArb8.h:103
TGeoTrap::fTl1
Double_t fTl1
Definition TGeoArb8.h:106
TGeoTrap::fH2
Double_t fH2
Definition TGeoArb8.h:108
TGeoTrap::Divide
TGeoVolume * Divide(TGeoVolume *voldiv, const char *divname, Int_t iaxis, Int_t ndiv, Double_t start, Double_t step) override
TGeoTrap::DistFromOutside_v
void DistFromOutside_v(const Double_t *points, const Double_t *dirs, Double_t *dists, Int_t vecsize, Double_t *step) const override
TGeoTrap::fTheta
Double_t fTheta
Definition TGeoArb8.h:102
TGeoTrap::GetMakeRuntimeShape
TGeoShape * GetMakeRuntimeShape(TGeoShape *mother, TGeoMatrix *mat) const override
TGeoTrap::DistFromOutside
Double_t DistFromOutside(const Double_t *point, const Double_t *dir, Int_t iact=1, Double_t step=TGeoShape::Big(), Double_t *safe=nullptr) const override
TGeoTrap::fAlpha1
Double_t fAlpha1
Definition TGeoArb8.h:107
TGeoVolumeMulti::AddVolume
void AddVolume(TGeoVolume *vol)
Add a volume with valid shape to the list of volumes.
Definition TGeoVolume.cxx:2709
TGeoVolume
TGeoVolume, TGeoVolumeMulti, TGeoVolumeAssembly are the volume classes.
Definition TGeoVolume.h:43
TGeoVolume::AddNodeOffset
void AddNodeOffset(TGeoVolume *vol, Int_t copy_no, Double_t offset=0, Option_t *option="")
Add a division node to the list of nodes.
Definition TGeoVolume.cxx:1138
TGeoVolume::GetMedium
TGeoMedium * GetMedium() const
Definition TGeoVolume.h:176
TGeoVolume::SetFinder
void SetFinder(TGeoPatternFinder *finder)
Definition TGeoVolume.h:245
TGeoVolume::GetNdaughters
Int_t GetNdaughters() const
Definition TGeoVolume.h:363
TGeoVolume::GetNodes
TObjArray * GetNodes()
Definition TGeoVolume.h:170
TObjArray::At
TObject * At(Int_t idx) const override
Definition TObjArray.h:170
TObject::TestBit
Bool_t TestBit(UInt_t f) const
Definition TObject.h:204
TObject::ClassName
virtual const char * ClassName() const
Returns name of class to which the object belongs.
Definition TObject.cxx:227
TObject::SetBit
void SetBit(UInt_t f, Bool_t set)
Set or unset the user status bits as specified in f.
Definition TObject.cxx:888
TObject::Error
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition TObject.cxx:1098
TObject::Fatal
virtual void Fatal(const char *method, const char *msgfmt,...) const
Issue fatal error message.
Definition TObject.cxx:1126
TString::Data
const char * Data() const
Definition TString.h:384
pt
TPaveText * pt
Definition entrylist_figure1.C:7
y
Double_t y[n]
Definition legend1.C:17
x
Double_t x[n]
Definition legend1.C:17
SetName
gr SetName("gr")
h1
TH1F * h1
Definition legend1.C:5
ROOT::ROOT_MATH_ARCH::gv_detail::dist
double dist(Rotation3D const &r1, Rotation3D const &r2)
Definition 3DDistances.cxx:42
TMath::LocMin
Long64_t LocMin(Long64_t n, const T *a)
Returns index of array with the minimum element.
Definition TMath.h:993
TMath::Max
Short_t Max(Short_t a, Short_t b)
Returns the largest of a and b.
Definition TMathBase.h:249
TMath::Sign
T1 Sign(T1 a, T2 b)
Returns a value with the magnitude of a and the sign of b.
Definition TMathBase.h:174
TMath::LocMax
Long64_t LocMax(Long64_t n, const T *a)
Returns index of array with the maximum element.
Definition TMath.h:1099
TMath::DegToRad
constexpr Double_t DegToRad()
Conversion from degree to radian: .
Definition TMath.h:82
TMath::Sqrt
Double_t Sqrt(Double_t x)
Returns the square root of x.
Definition TMath.h:673
TMath::Min
Short_t Min(Short_t a, Short_t b)
Returns the smallest of a and b.
Definition TMathBase.h:197
TMath::Cos
Double_t Cos(Double_t)
Returns the cosine of an angle of x radians.
Definition TMath.h:605
TMath::Sin
Double_t Sin(Double_t)
Returns the sine of an angle of x radians.
Definition TMath.h:599
TMath::Tan
Double_t Tan(Double_t)
Returns the tangent of an angle of x radians.
Definition TMath.h:611
TMath::Abs
Short_t Abs(Short_t d)
Returns the absolute value of parameter Short_t d.
Definition TMathBase.h:122
bvh::v2::cross
BVH_ALWAYS_INLINE Vec< T, 3 > cross(const Vec< T, 3 > &a, const Vec< T, 3 > &b)
Definition vec.h:104
hist015_TH1_read_and_draw_uhi.dz
dz
Definition hist015_TH1_read_and_draw_uhi.py:57
hist015_TH1_read_and_draw_uhi.ax
ax
Definition hist015_TH1_read_and_draw_uhi.py:40
hist015_TH1_read_and_draw_uhi.dx
dx
Definition hist015_TH1_read_and_draw_uhi.py:58
hist015_TH1_read_and_draw_uhi.dy
dy
Definition hist015_TH1_read_and_draw_uhi.py:59
v2
@ v2
Definition rootcling_impl.cxx:3557
v1
@ v1
Definition rootcling_impl.cxx:3556