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