Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
TGeoTorus.cxx
Go to the documentation of this file.
1// @(#)root/geom:$Id$
2// Author: Andrei Gheata 28/07/03
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/** \class TGeoTorus
13\ingroup Shapes_classes
14
15The torus is defined by its axial radius, its inner and outer radius.
16
17Begin_Macro
18{
19 TCanvas *c = new TCanvas("c", "c",0,0,600,600);
20 new TGeoManager("torus", "poza2");
21 TGeoMaterial *mat = new TGeoMaterial("Al", 26.98,13,2.7);
22 TGeoMedium *med = new TGeoMedium("MED",1,mat);
23 TGeoVolume *top = gGeoManager->MakeBox("TOP",med,100,100,100);
24 gGeoManager->SetTopVolume(top);
25 TGeoVolume *vol = gGeoManager->MakeTorus("TORUS",med, 40,20,25,0,270);
26 top->AddNode(vol,1);
27 gGeoManager->CloseGeometry();
28 gGeoManager->SetNsegments(30);
29 top->Draw();
30 TView *view = gPad->GetView();
31 if (view) view->ShowAxis();
32}
33End_Macro
34
35It may have a `phi `range:
36
37~~~{.cpp}
38TGeoTorus(Double_t R,Double_t Rmin,Double_t Rmax,Double_t Phi1,
39Double_t Dphi);
40~~~
41
42 - `R:` axial radius of the torus
43 - `Rmin:` inner radius
44 - `Rmax:` outer radius
45 - `Phi1:` starting phi angle
46 - `Dphi:` total phi range
47
48*/
49
50#include <iostream>
51
52#include "TGeoManager.h"
53#include "TGeoVolume.h"
54#include "TGeoTube.h"
55#include "TVirtualGeoPainter.h"
56#include "TGeoTorus.h"
57#include "TBuffer3D.h"
58#include "TBuffer3DTypes.h"
59#include "TMath.h"
60
61
62////////////////////////////////////////////////////////////////////////////////
63/// Default constructor
64
66{
68 fR = 0.0;
69 fRmin = 0.0;
70 fRmax = 0.0;
71 fPhi1 = 0.0;
72 fDphi = 0.0;
73}
74
75////////////////////////////////////////////////////////////////////////////////
76/// Constructor without name.
77
86
87////////////////////////////////////////////////////////////////////////////////
88/// Constructor with name.
89
99
100////////////////////////////////////////////////////////////////////////////////
101/// Constructor based on an array of parameters.
102/// - param[0] = R
103/// - param[1] = Rmin
104/// - param[2] = Rmax
105/// - param[3] = Phi1
106/// - param[4] = Dphi
107
109{
111 SetDimensions(param);
112 if (fRmin < 0 || fRmax < 0)
114 ComputeBBox();
115}
116
117////////////////////////////////////////////////////////////////////////////////
118/// Computes capacity of the shape in [length^3]
119
121{
122 Double_t capacity = (fDphi / 180.) * TMath::Pi() * TMath::Pi() * fR * (fRmax * fRmax - fRmin * fRmin);
123 return capacity;
124}
125
126////////////////////////////////////////////////////////////////////////////////
127/// Compute bounding box of the torus.
128
130{
131 fDZ = fRmax;
133 fDX = fDY = fR + fRmax;
134 return;
135 }
136 Double_t xc[4];
137 Double_t yc[4];
138 xc[0] = (fR + fRmax) * TMath::Cos(fPhi1 * TMath::DegToRad());
139 yc[0] = (fR + fRmax) * TMath::Sin(fPhi1 * TMath::DegToRad());
140 xc[1] = (fR + fRmax) * TMath::Cos((fPhi1 + fDphi) * TMath::DegToRad());
141 yc[1] = (fR + fRmax) * TMath::Sin((fPhi1 + fDphi) * TMath::DegToRad());
142 xc[2] = (fR - fRmax) * TMath::Cos(fPhi1 * TMath::DegToRad());
143 yc[2] = (fR - fRmax) * TMath::Sin(fPhi1 * TMath::DegToRad());
144 xc[3] = (fR - fRmax) * TMath::Cos((fPhi1 + fDphi) * TMath::DegToRad());
145 yc[3] = (fR - fRmax) * TMath::Sin((fPhi1 + fDphi) * TMath::DegToRad());
146
147 Double_t xmin = xc[TMath::LocMin(4, &xc[0])];
148 Double_t xmax = xc[TMath::LocMax(4, &xc[0])];
149 Double_t ymin = yc[TMath::LocMin(4, &yc[0])];
150 Double_t ymax = yc[TMath::LocMax(4, &yc[0])];
151 Double_t ddp = -fPhi1;
152 if (ddp < 0)
153 ddp += 360;
154 if (ddp <= fDphi)
155 xmax = fR + fRmax;
156 ddp = 90 - fPhi1;
157 if (ddp < 0)
158 ddp += 360;
159 if (ddp > 360)
160 ddp -= 360;
161 if (ddp <= fDphi)
162 ymax = fR + fRmax;
163 ddp = 180 - fPhi1;
164 if (ddp < 0)
165 ddp += 360;
166 if (ddp > 360)
167 ddp -= 360;
168 if (ddp <= fDphi)
169 xmin = -(fR + fRmax);
170 ddp = 270 - fPhi1;
171 if (ddp < 0)
172 ddp += 360;
173 if (ddp > 360)
174 ddp -= 360;
175 if (ddp <= fDphi)
176 ymin = -(fR + fRmax);
177 fOrigin[0] = (xmax + xmin) / 2;
178 fOrigin[1] = (ymax + ymin) / 2;
179 fOrigin[2] = 0;
180 fDX = (xmax - xmin) / 2;
181 fDY = (ymax - ymin) / 2;
182}
183
184////////////////////////////////////////////////////////////////////////////////
185/// Compute normal to closest surface from POINT.
186
187void TGeoTorus::ComputeNormal(const Double_t *point, const Double_t *dir, Double_t *norm) const
188{
189 Double_t phi = TMath::ATan2(point[1], point[0]);
190 if (fDphi < 360) {
197
198 Double_t daxis = Daxis(point, dir, 0);
199 if ((fRmax - daxis) > 1E-5) {
200 if (TGeoShape::IsSameWithinTolerance(fRmin, 0) || (daxis - fRmin) > 1E-5) {
201 TGeoShape::NormalPhi(point, dir, norm, c1, s1, c2, s2);
202 return;
203 }
204 }
205 }
206 Double_t r0[3];
207 r0[0] = fR * TMath::Cos(phi);
208 r0[1] = fR * TMath::Sin(phi);
209 r0[2] = 0;
210 Double_t normsq = 0;
211 for (Int_t i = 0; i < 3; i++) {
212 norm[i] = point[i] - r0[i];
213 normsq += norm[i] * norm[i];
214 }
215
217 norm[0] /= normsq;
218 norm[1] /= normsq;
219 norm[2] /= normsq;
220 if (dir[0] * norm[0] + dir[1] * norm[1] + dir[2] * norm[2] < 0) {
221 norm[0] = -norm[0];
222 norm[1] = -norm[1];
223 norm[2] = -norm[2];
224 }
225}
226
227////////////////////////////////////////////////////////////////////////////////
228/// Test if point is inside the torus.
229/// check phi range
230
232{
234 Double_t phi = TMath::ATan2(point[1], point[0]) * TMath::RadToDeg();
235 if (phi < 0)
236 phi += 360.0;
237 Double_t ddp = phi - fPhi1;
238 if (ddp < 0)
239 ddp += 360.;
240 if (ddp > fDphi)
241 return kFALSE;
242 }
243 // check radius
244 Double_t rxy = TMath::Sqrt(point[0] * point[0] + point[1] * point[1]);
245 Double_t radsq = (rxy - fR) * (rxy - fR) + point[2] * point[2];
246 if (radsq < fRmin * fRmin)
247 return kFALSE;
248 if (radsq > fRmax * fRmax)
249 return kFALSE;
250 return kTRUE;
251}
252
253////////////////////////////////////////////////////////////////////////////////
254/// Compute closest distance from point px,py to each vertex.
255
257{
259 Int_t numPoints = n * (n - 1);
260 if (fRmin > 0)
261 numPoints *= 2;
262 else if (fDphi < 360)
263 numPoints += 2;
264 return ShapeDistancetoPrimitive(numPoints, px, py);
265}
266
267////////////////////////////////////////////////////////////////////////////////
268/// Computes distance to axis of the torus from point pt + t*dir;
269
271{
272 Double_t p[3];
273 for (Int_t i = 0; i < 3; i++)
274 p[i] = pt[i] + t * dir[i];
275 Double_t rxy = TMath::Sqrt(p[0] * p[0] + p[1] * p[1]);
276 return TMath::Sqrt((rxy - fR) * (rxy - fR) + p[2] * p[2]);
277}
278
279////////////////////////////////////////////////////////////////////////////////
280/// Computes derivative w.r.t. t of the distance to axis of the torus from point pt + t*dir;
281
283{
284 Double_t p[3];
285 for (Int_t i = 0; i < 3; i++)
286 p[i] = pt[i] + t * dir[i];
287 Double_t rxy = TMath::Sqrt(p[0] * p[0] + p[1] * p[1]);
288 if (rxy < 1E-4)
289 return ((p[2] * dir[2] - fR * TMath::Sqrt(dir[0] * dir[0] + dir[1] * dir[1])) /
290 TMath::Sqrt(fR * fR + p[2] * p[2]));
291 Double_t d = TMath::Sqrt((rxy - fR) * (rxy - fR) + p[2] * p[2]);
293 return 0.;
294 Double_t dd = (p[0] * dir[0] + p[1] * dir[1] + p[2] * dir[2] - (p[0] * dir[0] + p[1] * dir[1]) * fR / rxy) / d;
295 return dd;
296}
297
298////////////////////////////////////////////////////////////////////////////////
299/// Second derivative of distance to torus axis w.r.t t.
300
302{
303 Double_t p[3];
304 for (Int_t i = 0; i < 3; i++)
305 p[i] = pt[i] + t * dir[i];
306 Double_t rxy = TMath::Sqrt(p[0] * p[0] + p[1] * p[1]);
307 if (rxy < 1E-6)
308 return 0;
309 Double_t daxis = TMath::Sqrt((rxy - fR) * (rxy - fR) + p[2] * p[2]);
311 return 0;
313 (p[0] * dir[0] + p[1] * dir[1] + p[2] * dir[2] - (p[0] * dir[0] + p[1] * dir[1]) * fR / rxy) / daxis;
314 Double_t dddaxis = 1 - ddaxis * ddaxis - (1 - dir[2] * dir[2]) * fR / rxy +
315 fR * (p[0] * dir[0] + p[1] * dir[1]) * (p[0] * dir[0] + p[1] * dir[1]) / (rxy * rxy * rxy);
316 dddaxis /= daxis;
317 return dddaxis;
318}
319
320////////////////////////////////////////////////////////////////////////////////
321/// Compute distance from inside point to surface of the torus.
322
325{
326 if (iact < 3 && safe) {
327 *safe = Safety(point, kTRUE);
328 if (iact == 0)
329 return TGeoShape::Big();
330 if ((iact == 1) && (step <= *safe))
331 return TGeoShape::Big();
332 }
333 Bool_t hasphi = (fDphi < 360);
334 Bool_t hasrmin = (fRmin > 0);
335 Double_t dout = ToBoundary(point, dir, fRmax, kTRUE);
336 // Double_t dax = Daxis(point,dir,dout);
337 Double_t din = (hasrmin) ? ToBoundary(point, dir, fRmin, kTRUE) : TGeoShape::Big();
339 if (snext > 1E10)
340 return TGeoShape::Tolerance();
341 if (hasphi) {
342 // Torus segment case.
343 Double_t c1, s1, c2, s2, cm, sm, cdfi;
346 c1 = TMath::Cos(phi1);
347 s1 = TMath::Sin(phi1);
348 c2 = TMath::Cos(phi2);
349 s2 = TMath::Sin(phi2);
350 Double_t fio = 0.5 * (phi1 + phi2);
351 cm = TMath::Cos(fio);
352 sm = TMath::Sin(fio);
353 cdfi = TMath::Cos(0.5 * (phi2 - phi1));
354 Double_t dphi =
355 TGeoTubeSeg::DistFromInsideS(point, dir, fR - fRmax, fR + fRmax, fRmax, c1, s1, c2, s2, cm, sm, cdfi);
356 Double_t daxis = Daxis(point, dir, dphi);
357 if (daxis >= fRmin + 1.E-8 && daxis <= fRmax - 1.E-8)
359 }
360 return snext;
361}
362
363////////////////////////////////////////////////////////////////////////////////
364/// Compute distance from outside point to surface of the torus.
365
368{
369 if (iact < 3 && safe) {
370 *safe = Safety(point, kFALSE);
371 if (iact == 0)
372 return TGeoShape::Big();
373 if ((iact == 1) && (step <= *safe))
374 return TGeoShape::Big();
375 }
376 // Check if the bounding box is crossed within the requested distance
377 Double_t sdist = TGeoBBox::DistFromOutside(point, dir, fDX, fDY, fDZ, fOrigin, step);
378 if (sdist >= step)
379 return TGeoShape::Big();
381 Bool_t hasphi = (fDphi < 360) ? kTRUE : kFALSE;
382 // Bool_t hasrmin = (fRmin>0)?kTRUE:kFALSE;
383 Double_t c1 = 0, s1 = 0, c2 = 0, s2 = 0, cm = 0, sm = 0, cdfi = 0;
385 Double_t phi, ddp, phi1, phi2, fio;
386 Double_t rxy2, dd;
388 Double_t pt[3];
389 Int_t i;
390
391 if (hasphi) {
392 // Torus segment case.
393 phi = TMath::ATan2(point[1], point[0]) * TMath::RadToDeg();
394 ;
395 if (phi < 0)
396 phi += 360;
397 ddp = phi - fPhi1;
398 if (ddp < 0)
399 ddp += 360;
400 ;
401 if (ddp <= fDphi)
402 inphi = kTRUE;
405 c1 = TMath::Cos(phi1);
406 s1 = TMath::Sin(phi1);
407 c2 = TMath::Cos(phi2);
408 s2 = TMath::Sin(phi2);
409 fio = 0.5 * (phi1 + phi2);
410 cm = TMath::Cos(fio);
411 sm = TMath::Sin(fio);
412 cdfi = TMath::Cos(0.5 * (phi2 - phi1));
413 }
414 // Check if we are inside or outside the bounding ring.
416 if (TMath::Abs(point[2]) <= fRmax) {
417 rxy2 = point[0] * point[0] + point[1] * point[1];
418 if ((rxy2 >= (fR - fRmax) * (fR - fRmax)) && (rxy2 <= (fR + fRmax) * (fR + fRmax))) {
419 if (!hasphi || inphi)
420 inbring = kTRUE;
421 }
422 }
423
424 // If outside the ring, compute distance to it.
426 Double_t eps = 1.E-8;
427 snext = 0;
428 daxis = -1;
429 memcpy(pt, point, 3 * sizeof(Double_t));
430 if (!inbring) {
431 if (hasphi)
432 dring = TGeoTubeSeg::DistFromOutsideS(point, dir, TMath::Max(0., fR - fRmax - eps), fR + fRmax + eps,
433 fRmax + eps, c1, s1, c2, s2, cm, sm, cdfi);
434 else
435 dring =
436 TGeoTube::DistFromOutsideS(point, dir, TMath::Max(0., fR - fRmax - eps), fR + fRmax + eps, fRmax + eps);
437 // If not crossing it, return BIG.
438 if (dring > 1E10)
439 return TGeoShape::Big();
440 snext = dring;
441 // Check if the crossing is due to phi.
442 daxis = Daxis(point, dir, snext);
443 if (daxis >= fRmin && daxis < fRmax)
444 return snext;
445 // Not a phi crossing -> propagate until we cross the ring.
446 for (i = 0; i < 3; i++)
447 pt[i] = point[i] + snext * dir[i];
448 }
449 // Point pt is inside the bounding ring, no phi crossing so far.
450 // Check if we are in the hole.
451 if (daxis < 0)
452 daxis = Daxis(pt, dir, 0);
453 if (daxis < fRmin + 1.E-8) {
454 // We are in the hole. Check if we came from outside.
455 if (snext > 0) {
456 // we can cross either the inner torus or exit the other hole.
457 snext += 0.1 * eps;
458 for (i = 0; i < 3; i++)
459 pt[i] += 0.1 * eps * dir[i];
460 }
461 // We are in the hole from the beginning.
462 // find first crossing with inner torus
463 dd = ToBoundary(pt, dir, fRmin, kFALSE);
464 // find exit distance from inner bounding ring
465 if (hasphi)
467 else
469 if (dd < dring)
470 return (snext + dd);
471 // we were exiting a hole inside phi hole
472 snext += dring + eps;
473 for (i = 0; i < 3; i++)
474 pt[i] = point[i] + snext * dir[i];
475 snext += DistFromOutside(pt, dir, 3);
476 return snext;
477 }
478 // We are inside the outer ring, having daxis>fRmax
479 // Compute distance to exit the bounding ring (again)
480 if (snext > 0) {
481 // we can cross either the inner torus or exit the other hole.
482 snext += 0.1 * eps;
483 for (i = 0; i < 3; i++)
484 pt[i] += 0.1 * eps * dir[i];
485 }
486 // Check intersection with outer torus
487 dd = ToBoundary(pt, dir, fRmax, kFALSE);
488 if (hasphi)
489 dring = TGeoTubeSeg::DistFromInsideS(pt, dir, TMath::Max(0., fR - fRmax - eps), fR + fRmax + eps, fRmax + eps, c1,
490 s1, c2, s2, cm, sm, cdfi);
491 else
492 dring = TGeoTube::DistFromInsideS(pt, dir, TMath::Max(0., fR - fRmax - eps), fR + fRmax + eps, fRmax + eps);
493 if (dd < dring) {
494 snext += dd;
495 return snext;
496 }
497 // We are exiting the bounding ring before crossing the torus -> propagate
498 snext += dring + eps;
499 for (i = 0; i < 3; i++)
500 pt[i] = point[i] + snext * dir[i];
501 snext += DistFromOutside(pt, dir, 3);
502 return snext;
503}
504
505////////////////////////////////////////////////////////////////////////////////
506/// Divide this torus shape belonging to volume "voldiv" into ndiv volumes
507/// called divname, from start position with the given step.
508
509TGeoVolume *TGeoTorus::Divide(TGeoVolume * /*voldiv*/, const char * /*divname*/, Int_t /*iaxis*/, Int_t /*ndiv*/,
510 Double_t /*start*/, Double_t /*step*/)
511{
512 return nullptr;
513}
514
515////////////////////////////////////////////////////////////////////////////////
516/// Returns name of axis IAXIS.
517
519{
520 switch (iaxis) {
521 case 1: return "R";
522 case 2: return "PHI";
523 case 3: return "Z";
524 default: return "UNDEFINED";
525 }
526}
527
528////////////////////////////////////////////////////////////////////////////////
529/// Get range of shape for a given axis.
530
532{
533 xlo = 0;
534 xhi = 0;
535 Double_t dx = 0;
536 switch (iaxis) {
537 case 1:
538 xlo = fRmin;
539 xhi = fRmax;
540 dx = xhi - xlo;
541 return dx;
542 case 2:
543 xlo = fPhi1;
544 xhi = fPhi1 + fDphi;
545 dx = fDphi;
546 return dx;
547 case 3: dx = 0; return dx;
548 }
549 return dx;
550}
551
552////////////////////////////////////////////////////////////////////////////////
553/// Fill vector param[4] with the bounding cylinder parameters. The order
554/// is the following : Rmin, Rmax, Phi1, Phi2, dZ
555
557{
558 param[0] = (fR - fRmax); // Rmin
559 param[1] = (fR + fRmax); // Rmax
560 param[2] = fPhi1; // Phi1
561 param[3] = fPhi1 + fDphi; // Phi2
562}
563
564////////////////////////////////////////////////////////////////////////////////
565/// Create a shape fitting the mother.
566
568{
570 return nullptr;
571 Error("GetMakeRuntimeShape", "parametrized toruses not supported");
572 return nullptr;
573}
574
575////////////////////////////////////////////////////////////////////////////////
576/// print shape parameters
577
579{
580 printf("*** Shape %s: TGeoTorus ***\n", GetName());
581 printf(" R = %11.5f\n", fR);
582 printf(" Rmin = %11.5f\n", fRmin);
583 printf(" Rmax = %11.5f\n", fRmax);
584 printf(" Phi1 = %11.5f\n", fPhi1);
585 printf(" Dphi = %11.5f\n", fDphi);
586 printf(" Bounding box:\n");
588}
589
590////////////////////////////////////////////////////////////////////////////////
591/// Creates a TBuffer3D describing *this* shape.
592/// Coordinates are in local reference frame.
593
595{
597 Int_t nbPnts = n * (n - 1);
598 Bool_t hasrmin = (GetRmin() > 0) ? kTRUE : kFALSE;
599 Bool_t hasphi = (GetDphi() < 360) ? kTRUE : kFALSE;
600 if (hasrmin)
601 nbPnts *= 2;
602 else if (hasphi)
603 nbPnts += 2;
604
605 Int_t nbSegs = (2 * n - 1) * (n - 1);
606 Int_t nbPols = (n - 1) * (n - 1);
607 if (hasrmin) {
608 nbSegs += (2 * n - 1) * (n - 1);
609 nbPols += (n - 1) * (n - 1);
610 }
611 if (hasphi) {
612 nbSegs += 2 * (n - 1);
613 nbPols += 2 * (n - 1);
614 }
615
616 TBuffer3D *buff =
618 if (buff) {
619 SetPoints(buff->fPnts);
621 }
622
623 return buff;
624}
625
626////////////////////////////////////////////////////////////////////////////////
627/// Fill TBuffer3D structure for segments and polygons.
628
630{
631 Int_t i, j;
633 // Int_t nbPnts = n*(n-1);
634 Int_t indx, indp, startcap = 0;
635 Bool_t hasrmin = (GetRmin() > 0) ? kTRUE : kFALSE;
636 Bool_t hasphi = (GetDphi() < 360) ? kTRUE : kFALSE;
637 // if (hasrmin) nbPnts *= 2;
638 // else if (hasphi) nbPnts += 2;
640
641 indp = n * (n - 1); // start index for points on inner surface
642 memset(buff.fSegs, 0, buff.NbSegs() * 3 * sizeof(Int_t));
643
644 // outer surface phi circles = n*(n-1) -> [0, n*(n-1) -1]
645 // connect point j with point j+1 on same row
646 indx = 0;
647 for (i = 0; i < n; i++) { // rows [0,n-1]
648 for (j = 0; j < n - 1; j++) { // points on a row [0, n-2]
649 buff.fSegs[indx + (i * (n - 1) + j) * 3] = c;
650 buff.fSegs[indx + (i * (n - 1) + j) * 3 + 1] = i * (n - 1) + j; // j on row i
651 buff.fSegs[indx + (i * (n - 1) + j) * 3 + 2] = i * (n - 1) + ((j + 1) % (n - 1)); // j+1 on row i
652 }
653 }
654 indx += 3 * n * (n - 1);
655 // outer surface generators = (n-1)*(n-1) -> [n*(n-1), (2*n-1)*(n-1) -1]
656 // connect point j on row i with point j on row i+1
657 for (i = 0; i < n - 1; i++) { // rows [0, n-2]
658 for (j = 0; j < n - 1; j++) { // points on a row [0, n-2]
659 buff.fSegs[indx + (i * (n - 1) + j) * 3] = c;
660 buff.fSegs[indx + (i * (n - 1) + j) * 3 + 1] = i * (n - 1) + j; // j on row i
661 buff.fSegs[indx + (i * (n - 1) + j) * 3 + 2] = (i + 1) * (n - 1) + j; // j on row i+1
662 }
663 }
664 indx += 3 * (n - 1) * (n - 1);
665 startcap = (2 * n - 1) * (n - 1);
666
667 if (hasrmin) {
668 // inner surface phi circles = n*(n-1) -> [(2*n-1)*(n-1), (3*n-1)*(n-1) -1]
669 // connect point j with point j+1 on same row
670 for (i = 0; i < n; i++) { // rows [0, n-1]
671 for (j = 0; j < n - 1; j++) { // points on a row [0, n-2]
672 buff.fSegs[indx + (i * (n - 1) + j) * 3] = c; // lighter color
673 buff.fSegs[indx + (i * (n - 1) + j) * 3 + 1] = indp + i * (n - 1) + j; // j on row i
674 buff.fSegs[indx + (i * (n - 1) + j) * 3 + 2] = indp + i * (n - 1) + ((j + 1) % (n - 1)); // j+1 on row i
675 }
676 }
677 indx += 3 * n * (n - 1);
678 // inner surface generators = (n-1)*n -> [(3*n-1)*(n-1), (4*n-2)*(n-1) -1]
679 // connect point j on row i with point j on row i+1
680 for (i = 0; i < n - 1; i++) { // rows [0, n-2]
681 for (j = 0; j < n - 1; j++) { // points on a row [0, n-2]
682 buff.fSegs[indx + (i * (n - 1) + j) * 3] = c; // lighter color
683 buff.fSegs[indx + (i * (n - 1) + j) * 3 + 1] = indp + i * (n - 1) + j; // j on row i
684 buff.fSegs[indx + (i * (n - 1) + j) * 3 + 2] = indp + (i + 1) * (n - 1) + j; // j on row i+1
685 }
686 }
687 indx += 3 * (n - 1) * (n - 1);
688 startcap = (4 * n - 2) * (n - 1);
689 }
690
691 if (hasphi) {
692 if (hasrmin) {
693 // endcaps = 2*(n-1) -> [(4*n-2)*(n-1), 4*n*(n-1)-1]
694 i = 0;
695 for (j = 0; j < n - 1; j++) {
696 buff.fSegs[indx + j * 3] = c + 1;
697 buff.fSegs[indx + j * 3 + 1] = (n - 1) * i + j; // outer j on row 0
698 buff.fSegs[indx + j * 3 + 2] = indp + (n - 1) * i + j; // inner j on row 0
699 }
700 indx += 3 * (n - 1);
701 i = n - 1;
702 for (j = 0; j < n - 1; j++) {
703 buff.fSegs[indx + j * 3] = c + 1;
704 buff.fSegs[indx + j * 3 + 1] = (n - 1) * i + j; // outer j on row n-1
705 buff.fSegs[indx + j * 3 + 2] = indp + (n - 1) * i + j; // inner j on row n-1
706 }
707 indx += 3 * (n - 1);
708 } else {
709 i = 0;
710 for (j = 0; j < n - 1; j++) {
711 buff.fSegs[indx + j * 3] = c + 1;
712 buff.fSegs[indx + j * 3 + 1] = (n - 1) * i + j; // outer j on row 0
713 buff.fSegs[indx + j * 3 + 2] = n * (n - 1); // center of first endcap
714 }
715 indx += 3 * (n - 1);
716 i = n - 1;
717 for (j = 0; j < n - 1; j++) {
718 buff.fSegs[indx + j * 3] = c + 1;
719 buff.fSegs[indx + j * 3 + 1] = (n - 1) * i + j; // outer j on row n-1
720 buff.fSegs[indx + j * 3 + 2] = n * (n - 1) + 1; // center of second endcap
721 }
722 indx += 3 * (n - 1);
723 }
724 }
725
726 indx = 0;
727 memset(buff.fPols, 0, buff.NbPols() * 6 * sizeof(Int_t));
728
729 // outer surface = (n-1)*(n-1) -> [0, (n-1)*(n-1)-1]
730 // normal pointing out
731 for (i = 0; i < n - 1; i++) {
732 for (j = 0; j < n - 1; j++) {
733 buff.fPols[indx++] = c;
734 buff.fPols[indx++] = 4;
735 buff.fPols[indx++] = n * (n - 1) + (n - 1) * i + ((j + 1) % (n - 1)); // generator j+1 on outer row i
736 buff.fPols[indx++] = (n - 1) * (i + 1) + j; // seg j on outer row i+1
737 buff.fPols[indx++] = n * (n - 1) + (n - 1) * i + j; // generator j on outer row i
738 buff.fPols[indx++] = (n - 1) * i + j; // seg j on outer row i
739 }
740 }
741 if (hasrmin) {
742 indp = (2 * n - 1) * (n - 1); // start index of inner segments
743 // inner surface = (n-1)*(n-1) -> [(n-1)*(n-1), 2*(n-1)*(n-1)-1]
744 // normal pointing out
745 for (i = 0; i < n - 1; i++) {
746 for (j = 0; j < n - 1; j++) {
747 buff.fPols[indx++] = c;
748 buff.fPols[indx++] = 4;
749 buff.fPols[indx++] = indp + n * (n - 1) + (n - 1) * i + j; // generator j on inner row i
750 buff.fPols[indx++] = indp + (n - 1) * (i + 1) + j; // seg j on inner row i+1
751 buff.fPols[indx++] = indp + n * (n - 1) + (n - 1) * i + ((j + 1) % (n - 1)); // generator j+1 on inner r>
752 buff.fPols[indx++] = indp + (n - 1) * i + j; // seg j on inner row i
753 }
754 }
755 }
756 if (hasphi) {
757 // endcaps = 2*(n-1) -> [2*(n-1)*(n-1), 2*n*(n-1)-1]
758 i = 0; // row 0
759 Int_t np = (hasrmin) ? 4 : 3;
760 for (j = 0; j < n - 1; j++) {
761 buff.fPols[indx++] = c + 1;
762 buff.fPols[indx++] = np;
763 buff.fPols[indx++] = j; // seg j on outer row 0 a
764 buff.fPols[indx++] = startcap + j; // endcap j on row 0 d
765 if (hasrmin)
766 buff.fPols[indx++] = indp + j; // seg j on inner row 0 c
767 buff.fPols[indx++] = startcap + ((j + 1) % (n - 1)); // endcap j+1 on row 0 b
768 }
769
770 i = n - 1; // row n-1
771 for (j = 0; j < n - 1; j++) {
772 buff.fPols[indx++] = c + 1;
773 buff.fPols[indx++] = np;
774 buff.fPols[indx++] = (n - 1) * i + j; // seg j on outer row n-1 a
775 buff.fPols[indx++] = startcap + (n - 1) + ((j + 1) % (n - 1)); // endcap j+1 on row n-1 d
776 if (hasrmin)
777 buff.fPols[indx++] = indp + (n - 1) * i + j; // seg j on inner row n-1 c
778 buff.fPols[indx++] = startcap + (n - 1) + j; // endcap j on row n-1 b
779 }
780 }
781}
782
783////////////////////////////////////////////////////////////////////////////////
784/// computes the closest distance from given point to this shape, according
785/// to option. The matching point on the shape is stored in spoint.
786
788{
789 Double_t saf[2];
790 Int_t i;
791 Double_t rxy = TMath::Sqrt(point[0] * point[0] + point[1] * point[1]);
792 Double_t rad = TMath::Sqrt((rxy - fR) * (rxy - fR) + point[2] * point[2]);
793 saf[0] = rad - fRmin;
794 saf[1] = fRmax - rad;
796 if (in)
797 return TMath::Min(saf[0], saf[1]);
798 for (i = 0; i < 2; i++)
799 saf[i] = -saf[i];
800 return TMath::Max(saf[0], saf[1]);
801 }
802
804 if (in) {
805 Double_t safe = TMath::Min(saf[0], saf[1]);
806 return TMath::Min(safe, safphi);
807 }
808 for (i = 0; i < 2; i++)
809 saf[i] = -saf[i];
810 Double_t safe = TMath::Max(saf[0], saf[1]);
811 return TMath::Max(safe, safphi);
812}
813
814////////////////////////////////////////////////////////////////////////////////
815/// Save a primitive as a C++ statement(s) on output stream "out".
816
817void TGeoTorus::SavePrimitive(std::ostream &out, Option_t * /*option*/ /*= ""*/)
818{
820 return;
821 out << " // Shape: " << GetName() << " type: " << ClassName() << std::endl;
822 out << " r = " << fR << ";" << std::endl;
823 out << " rmin = " << fRmin << ";" << std::endl;
824 out << " rmax = " << fRmax << ";" << std::endl;
825 out << " phi1 = " << fPhi1 << ";" << std::endl;
826 out << " dphi = " << fDphi << ";" << std::endl;
827 out << " TGeoShape *" << GetPointerName() << " = new TGeoTorus(\"" << GetName() << "\",r,rmin,rmax,phi1,dphi);"
828 << std::endl;
830}
831
832////////////////////////////////////////////////////////////////////////////////
833/// Set torus dimensions.
834
836{
837 fR = r;
838 fRmin = rmin;
839 fRmax = rmax;
840 fPhi1 = phi1;
841 if (fPhi1 < 0)
842 fPhi1 += 360.;
843 fDphi = dphi;
844}
845
846////////////////////////////////////////////////////////////////////////////////
847/// Set torus dimensions starting from a list.
848
850{
851 SetTorusDimensions(param[0], param[1], param[2], param[3], param[4]);
852}
853
854////////////////////////////////////////////////////////////////////////////////
855/// Create torus mesh points
856
858{
859 if (!points)
860 return;
863 Double_t dpin = 360. / (n - 1);
864 Double_t dpout = fDphi / (n - 1);
865 Double_t co, so, ci, si;
867 Int_t i, j;
868 Int_t indx = 0;
869 // loop outer mesh -> n*(n-1) points [0, 3*n*(n-1)-1]
870 for (i = 0; i < n; i++) {
871 phout = (fPhi1 + i * dpout) * TMath::DegToRad();
873 so = TMath::Sin(phout);
874 for (j = 0; j < n - 1; j++) {
875 phin = j * dpin * TMath::DegToRad();
876 ci = TMath::Cos(phin);
877 si = TMath::Sin(phin);
878 points[indx++] = (fR + fRmax * ci) * co;
879 points[indx++] = (fR + fRmax * ci) * so;
880 points[indx++] = fRmax * si;
881 }
882 }
883
884 if (havermin) {
885 // loop inner mesh -> n*(n-1) points [3*n*(n-1), 6*n*(n-1)]
886 for (i = 0; i < n; i++) {
887 phout = (fPhi1 + i * dpout) * TMath::DegToRad();
889 so = TMath::Sin(phout);
890 for (j = 0; j < n - 1; j++) {
891 phin = j * dpin * TMath::DegToRad();
892 ci = TMath::Cos(phin);
893 si = TMath::Sin(phin);
894 points[indx++] = (fR + fRmin * ci) * co;
895 points[indx++] = (fR + fRmin * ci) * so;
896 points[indx++] = fRmin * si;
897 }
898 }
899 } else {
900 if (fDphi < 360.) {
901 // just add extra 2 points on the centers of the 2 phi cuts [3*n*n, 3*n*n+1]
904 points[indx++] = 0;
907 points[indx++] = 0;
908 }
909 }
910}
911
912////////////////////////////////////////////////////////////////////////////////
913/// Create torus mesh points
914
916{
917 if (!points)
918 return;
921 Double_t dpin = 360. / (n - 1);
922 Double_t dpout = fDphi / (n - 1);
923 Double_t co, so, ci, si;
925 Int_t i, j;
926 Int_t indx = 0;
927 // loop outer mesh -> n*(n-1) points [0, 3*n*(n-1)-1]
928 // plane i = 0, n-1 point j = 0, n-1 ipoint = n*i + j
929 for (i = 0; i < n; i++) {
930 phout = (fPhi1 + i * dpout) * TMath::DegToRad();
932 so = TMath::Sin(phout);
933 for (j = 0; j < n - 1; j++) {
934 phin = j * dpin * TMath::DegToRad();
935 ci = TMath::Cos(phin);
936 si = TMath::Sin(phin);
937 points[indx++] = (fR + fRmax * ci) * co;
938 points[indx++] = (fR + fRmax * ci) * so;
939 points[indx++] = fRmax * si;
940 }
941 }
942
943 if (havermin) {
944 // loop inner mesh -> n*(n-1) points [3*n*(n-1), 6*n*(n-1)]
945 // plane i = 0, n-1 point j = 0, n-1 ipoint = n*n + n*i + j
946 for (i = 0; i < n; i++) {
947 phout = (fPhi1 + i * dpout) * TMath::DegToRad();
949 so = TMath::Sin(phout);
950 for (j = 0; j < n - 1; j++) {
951 phin = j * dpin * TMath::DegToRad();
952 ci = TMath::Cos(phin);
953 si = TMath::Sin(phin);
954 points[indx++] = (fR + fRmin * ci) * co;
955 points[indx++] = (fR + fRmin * ci) * so;
956 points[indx++] = fRmin * si;
957 }
958 }
959 } else {
960 if (fDphi < 360.) {
961 // just add extra 2 points on the centers of the 2 phi cuts [n*n, n*n+1]
962 // ip1 = n*(n-1) + 0;
963 // ip2 = n*(n-1) + 1
966 points[indx++] = 0;
969 points[indx++] = 0;
970 }
971 }
972}
973
974////////////////////////////////////////////////////////////////////////////////
975/// Return number of vertices of the mesh representation
976
978{
980 Int_t numPoints = n * (n - 1);
982 numPoints *= 2;
983 else if (fDphi < 360.)
984 numPoints += 2;
985 return numPoints;
986}
987
988////////////////////////////////////////////////////////////////////////////////
989/// fill size of this 3-D object
990
991void TGeoTorus::Sizeof3D() const {}
992
993////////////////////////////////////////////////////////////////////////////////
994/// Find real solutions of the cubic equation : x^3 + a*x^2 + b*x + c = 0
995/// Input: a,b,c
996/// Output: x[3] real solutions
997/// Returns number of real solutions (1 or 3)
998
1000{
1001 const Double_t ott = 1. / 3.;
1002 const Double_t sq3 = TMath::Sqrt(3.);
1003 Int_t ireal = 1;
1004 Double_t p = b - a * a * ott;
1005 Double_t q = c - a * b * ott + 2. * a * a * a * ott * ott * ott;
1006 Double_t delta = 4 * p * p * p + 27 * q * q;
1007 // Double_t y1r, y1i, y2r, y2i;
1008 Double_t t, u;
1009 if (delta >= 0) {
1010 delta = TMath::Sqrt(delta);
1011 t = (-3 * q * sq3 + delta) / (6 * sq3);
1012 u = (3 * q * sq3 + delta) / (6 * sq3);
1013 x[0] = TMath::Sign(1., t) * TMath::Power(TMath::Abs(t), ott) -
1014 TMath::Sign(1., u) * TMath::Power(TMath::Abs(u), ott) - a * ott;
1015 } else {
1016 delta = TMath::Sqrt(-delta);
1017 t = -0.5 * q;
1018 u = delta / (6 * sq3);
1019 x[0] = 2. * TMath::Power(t * t + u * u, 0.5 * ott) * TMath::Cos(ott * TMath::ATan2(u, t));
1020 x[0] -= a * ott;
1021 }
1022
1023 t = x[0] * x[0] + a * x[0] + b;
1024 u = a + x[0];
1025 delta = u * u - 4. * t;
1026 if (delta >= 0) {
1027 ireal = 3;
1028 delta = TMath::Sqrt(delta);
1029 x[1] = 0.5 * (-u - delta);
1030 x[2] = 0.5 * (-u + delta);
1031 }
1032 return ireal;
1033}
1034
1035////////////////////////////////////////////////////////////////////////////////
1036/// Find real solutions of the quartic equation : x^4 + a*x^3 + b*x^2 + c*x + d = 0
1037/// Input: a,b,c,d
1038/// Output: x[4] - real solutions
1039/// Returns number of real solutions (0 to 3)
1040
1042{
1043 Double_t e = b - 3. * a * a / 8.;
1044 Double_t f = c + a * a * a / 8. - 0.5 * a * b;
1045 Double_t g = d - 3. * a * a * a * a / 256. + a * a * b / 16. - a * c / 4.;
1046 Double_t xx[4];
1047 Int_t ind[4];
1048 Double_t delta;
1049 Double_t h = 0;
1050 Int_t ireal = 0;
1051 Int_t i;
1053 delta = e * e - 4. * g;
1054 if (delta < 0)
1055 return 0;
1056 delta = TMath::Sqrt(delta);
1057 h = 0.5 * (-e - delta);
1058 if (h >= 0) {
1059 h = TMath::Sqrt(h);
1060 x[ireal++] = -h - 0.25 * a;
1061 x[ireal++] = h - 0.25 * a;
1062 }
1063 h = 0.5 * (-e + delta);
1064 if (h >= 0) {
1065 h = TMath::Sqrt(h);
1066 x[ireal++] = -h - 0.25 * a;
1067 x[ireal++] = h - 0.25 * a;
1068 }
1069 if (ireal > 0) {
1071 for (i = 0; i < ireal; i++)
1072 xx[i] = x[ind[i]];
1073 memcpy(x, xx, ireal * sizeof(Double_t));
1074 }
1075 return ireal;
1076 }
1077
1079 x[ireal++] = -0.25 * a;
1080 ind[0] = SolveCubic(0, e, f, xx);
1081 for (i = 0; i < ind[0]; i++)
1082 x[ireal++] = xx[i] - 0.25 * a;
1083 if (ireal > 0) {
1085 for (i = 0; i < ireal; i++)
1086 xx[i] = x[ind[i]];
1087 memcpy(x, xx, ireal * sizeof(Double_t));
1088 }
1089 return ireal;
1090 }
1091
1092 ireal = SolveCubic(2. * e, e * e - 4. * g, -f * f, xx);
1093 if (ireal == 1) {
1094 if (xx[0] <= 0)
1095 return 0;
1096 h = TMath::Sqrt(xx[0]);
1097 } else {
1098 // 3 real solutions of the cubic
1099 for (i = 0; i < 3; i++) {
1100 h = xx[i];
1101 if (h >= 0)
1102 break;
1103 }
1104 if (h <= 0)
1105 return 0;
1106 h = TMath::Sqrt(h);
1107 }
1108 Double_t j = 0.5 * (e + h * h - f / h);
1109 ireal = 0;
1110 delta = h * h - 4. * j;
1111 if (delta >= 0) {
1112 delta = TMath::Sqrt(delta);
1113 x[ireal++] = 0.5 * (-h - delta) - 0.25 * a;
1114 x[ireal++] = 0.5 * (-h + delta) - 0.25 * a;
1115 }
1116 delta = h * h - 4. * g / j;
1117 if (delta >= 0) {
1118 delta = TMath::Sqrt(delta);
1119 x[ireal++] = 0.5 * (h - delta) - 0.25 * a;
1120 x[ireal++] = 0.5 * (h + delta) - 0.25 * a;
1121 }
1122 if (ireal > 0) {
1124 for (i = 0; i < ireal; i++)
1125 xx[i] = x[ind[i]];
1126 memcpy(x, xx, ireal * sizeof(Double_t));
1127 }
1128 return ireal;
1129}
1130
1131////////////////////////////////////////////////////////////////////////////////
1132/// Returns distance to the surface or the torus (fR,r) from a point, along
1133/// a direction. Point is close enough to the boundary so that the distance
1134/// to the torus is decreasing while moving along the given direction.
1135
1137{
1138 // Compute coefficients of the quartic
1140 Double_t r0sq = pt[0] * pt[0] + pt[1] * pt[1] + pt[2] * pt[2];
1141 Double_t rdotn = pt[0] * dir[0] + pt[1] * dir[1] + pt[2] * dir[2];
1142 Double_t rsumsq = fR * fR + r * r;
1143 Double_t a = 4. * rdotn;
1144 Double_t b = 2. * (r0sq + 2. * rdotn * rdotn - rsumsq + 2. * fR * fR * dir[2] * dir[2]);
1145 Double_t c = 4. * (r0sq * rdotn - rsumsq * rdotn + 2. * fR * fR * pt[2] * dir[2]);
1146 Double_t d = r0sq * r0sq - 2. * r0sq * rsumsq + 4. * fR * fR * pt[2] * pt[2] + (fR * fR - r * r) * (fR * fR - r * r);
1147
1148 Double_t x[4], y[4];
1149 Int_t nsol = 0;
1150
1151 if (TMath::Abs(dir[2]) < 1E-3 && TMath::Abs(pt[2]) < 0.1 * r) {
1152 Double_t r0 = fR - TMath::Sqrt((r - pt[2]) * (r + pt[2]));
1153 Double_t b0 = (pt[0] * dir[0] + pt[1] * dir[1]) / (dir[0] * dir[0] + dir[1] * dir[1]);
1154 Double_t c0 = (pt[0] * pt[0] + (pt[1] - r0) * (pt[1] + r0)) / (dir[0] * dir[0] + dir[1] * dir[1]);
1155 Double_t delta = b0 * b0 - c0;
1156 if (delta > 0) {
1157 y[nsol] = -b0 - TMath::Sqrt(delta);
1158 if (y[nsol] > -tol)
1159 nsol++;
1160 y[nsol] = -b0 + TMath::Sqrt(delta);
1161 if (y[nsol] > -tol)
1162 nsol++;
1163 }
1164 r0 = fR + TMath::Sqrt((r - pt[2]) * (r + pt[2]));
1165 c0 = (pt[0] * pt[0] + (pt[1] - r0) * (pt[1] + r0)) / (dir[0] * dir[0] + dir[1] * dir[1]);
1166 delta = b0 * b0 - c0;
1167 if (delta > 0) {
1168 y[nsol] = -b0 - TMath::Sqrt(delta);
1169 if (y[nsol] > -tol)
1170 nsol++;
1171 y[nsol] = -b0 + TMath::Sqrt(delta);
1172 if (y[nsol] > -tol)
1173 nsol++;
1174 }
1175 if (nsol) {
1176 // Sort solutions
1177 Int_t ind[4];
1179 for (Int_t j = 0; j < nsol; j++)
1180 x[j] = y[ind[j]];
1181 }
1182 } else {
1183 nsol = SolveQuartic(a, b, c, d, x);
1184 }
1185 if (!nsol)
1186 return TGeoShape::Big();
1187 // look for first positive solution
1188 Double_t phi, ndotd;
1189 Double_t r0[3], norm[3];
1191 for (Int_t i = 0; i < nsol; i++) {
1192 if (x[i] < -10)
1193 continue;
1194 phi = TMath::ATan2(pt[1] + x[i] * dir[1], pt[0] + x[i] * dir[0]);
1195 r0[0] = fR * TMath::Cos(phi);
1196 r0[1] = fR * TMath::Sin(phi);
1197 r0[2] = 0;
1198 for (Int_t ipt = 0; ipt < 3; ipt++)
1199 norm[ipt] = pt[ipt] + x[i] * dir[ipt] - r0[ipt];
1200 ndotd = norm[0] * dir[0] + norm[1] * dir[1] + norm[2] * dir[2];
1201 if (inner ^ in) {
1202 if (ndotd < 0)
1203 continue;
1204 } else {
1205 if (ndotd > 0)
1206 continue;
1207 }
1208 Double_t s = x[i];
1209 Double_t eps = TGeoShape::Big();
1210 Double_t delta = s * s * s * s + a * s * s * s + b * s * s + c * s + d;
1211 Double_t eps0 = -delta / (4. * s * s * s + 3. * a * s * s + 2. * b * s + c);
1212 while (TMath::Abs(eps) > TGeoShape::Tolerance()) {
1213 if (TMath::Abs(eps0) > 100)
1214 break;
1215 s += eps0;
1217 break;
1218 delta = s * s * s * s + a * s * s * s + b * s * s + c * s + d;
1219 eps = -delta / (4. * s * s * s + 3. * a * s * s + 2. * b * s + c);
1220 if (TMath::Abs(eps) > TMath::Abs(eps0))
1221 break;
1222 eps0 = eps;
1223 }
1224 if (s < -TGeoShape::Tolerance())
1225 continue;
1226 return TMath::Max(0., s);
1227 }
1228 return TGeoShape::Big();
1229}
1230
1231////////////////////////////////////////////////////////////////////////////////
1232/// Returns numbers of vertices, segments and polygons composing the shape mesh.
1233
1235{
1237 nvert = n * (n - 1);
1238 Bool_t hasrmin = (GetRmin() > 0) ? kTRUE : kFALSE;
1239 Bool_t hasphi = (GetDphi() < 360) ? kTRUE : kFALSE;
1240 if (hasrmin)
1241 nvert *= 2;
1242 else if (hasphi)
1243 nvert += 2;
1244 nsegs = (2 * n - 1) * (n - 1);
1245 npols = (n - 1) * (n - 1);
1246 if (hasrmin) {
1247 nsegs += (2 * n - 1) * (n - 1);
1248 npols += (n - 1) * (n - 1);
1249 }
1250 if (hasphi) {
1251 nsegs += 2 * (n - 1);
1252 npols += 2 * (n - 1);
1253 }
1254}
1255
1256////////////////////////////////////////////////////////////////////////////////
1257/// Fills a static 3D buffer and returns a reference.
1258
1260{
1261 static TBuffer3D buffer(TBuffer3DTypes::kGeneric);
1262
1264
1267 Int_t nbPnts = n * (n - 1);
1268 Bool_t hasrmin = (GetRmin() > 0) ? kTRUE : kFALSE;
1269 Bool_t hasphi = (GetDphi() < 360) ? kTRUE : kFALSE;
1270 if (hasrmin)
1271 nbPnts *= 2;
1272 else if (hasphi)
1273 nbPnts += 2;
1274
1275 Int_t nbSegs = (2 * n - 1) * (n - 1);
1276 Int_t nbPols = (n - 1) * (n - 1);
1277 if (hasrmin) {
1278 nbSegs += (2 * n - 1) * (n - 1);
1279 nbPols += (n - 1) * (n - 1);
1280 }
1281 if (hasphi) {
1282 nbSegs += 2 * (n - 1);
1283 nbPols += 2 * (n - 1);
1284 }
1285
1286 if (buffer.SetRawSizes(nbPnts, 3 * nbPnts, nbSegs, 3 * nbSegs, nbPols, 6 * nbPols)) {
1288 }
1289 }
1290 // TODO: Push down to TGeoShape?? But would have to do raw sizes set first..
1291 // can rest of TGeoShape be deferred until after
1293 SetPoints(buffer.fPnts);
1294 if (!buffer.fLocalFrame) {
1295 TransformPoints(buffer.fPnts, buffer.NbPnts());
1296 }
1297
1298 SetSegsAndPols(buffer);
1300 }
1301
1302 return buffer;
1303}
1304
1305////////////////////////////////////////////////////////////////////////////////
1306/// Check the inside status for each of the points in the array.
1307/// Input: Array of point coordinates + vector size
1308/// Output: Array of Booleans for the inside of each point
1309
1311{
1312 for (Int_t i = 0; i < vecsize; i++)
1313 inside[i] = Contains(&points[3 * i]);
1314}
1315
1316////////////////////////////////////////////////////////////////////////////////
1317/// Compute the normal for an array o points so that norm.dot.dir is positive
1318/// Input: Arrays of point coordinates and directions + vector size
1319/// Output: Array of normal directions
1320
1322{
1323 for (Int_t i = 0; i < vecsize; i++)
1324 ComputeNormal(&points[3 * i], &dirs[3 * i], &norms[3 * i]);
1325}
1326
1327////////////////////////////////////////////////////////////////////////////////
1328/// Compute distance from array of input points having directions specified by dirs. Store output in dists
1329
1331 Double_t *step) const
1332{
1333 for (Int_t i = 0; i < vecsize; i++)
1334 dists[i] = DistFromInside(&points[3 * i], &dirs[3 * i], 3, step[i]);
1335}
1336
1337////////////////////////////////////////////////////////////////////////////////
1338/// Compute distance from array of input points having directions specified by dirs. Store output in dists
1339
1341 Double_t *step) const
1342{
1343 for (Int_t i = 0; i < vecsize; i++)
1344 dists[i] = DistFromOutside(&points[3 * i], &dirs[3 * i], 3, step[i]);
1345}
1346
1347////////////////////////////////////////////////////////////////////////////////
1348/// Compute safe distance from each of the points in the input array.
1349/// Input: Array of point coordinates, array of statuses for these points, size of the arrays
1350/// Output: Safety values
1351
1353{
1354 for (Int_t i = 0; i < vecsize; i++)
1355 safe[i] = Safety(&points[3 * i], inside[i]);
1356}
#define d(i)
Definition RSha256.hxx:102
#define b(i)
Definition RSha256.hxx:100
#define f(i)
Definition RSha256.hxx:104
#define c(i)
Definition RSha256.hxx:101
#define g(i)
Definition RSha256.hxx:105
#define a(i)
Definition RSha256.hxx:99
#define s1(x)
Definition RSha256.hxx:91
#define h(i)
Definition RSha256.hxx:106
#define e(i)
Definition RSha256.hxx:103
float Float_t
Float 4 bytes (float)
Definition RtypesCore.h:71
constexpr Bool_t kFALSE
Definition RtypesCore.h:108
double Double_t
Double 8 bytes.
Definition RtypesCore.h:73
constexpr Bool_t kTRUE
Definition RtypesCore.h:107
const char Option_t
Option string (const char)
Definition RtypesCore.h:80
ROOT::Detail::TRangeCast< T, true > TRangeDynCast
TRangeDynCast is an adapter class that allows the typed iteration through a TCollection.
winID h TVirtualViewer3D TVirtualGLPainter p
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t np
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t r
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t points
char name[80]
Definition TGX11.cxx:110
R__EXTERN TGeoManager * gGeoManager
float xmin
float * q
float ymin
float xmax
float ymax
Generic 3D primitive description class.
Definition TBuffer3D.h:18
UInt_t NbPnts() const
Definition TBuffer3D.h:80
Bool_t SectionsValid(UInt_t mask) const
Definition TBuffer3D.h:67
void SetSectionsValid(UInt_t mask)
Definition TBuffer3D.h:65
Bool_t fLocalFrame
Definition TBuffer3D.h:90
Bool_t SetRawSizes(UInt_t reqPnts, UInt_t reqPntsCapacity, UInt_t reqSegs, UInt_t reqSegsCapacity, UInt_t reqPols, UInt_t reqPolsCapacity)
Set kRaw tessellation section of buffer with supplied sizes.
Double_t * fPnts
Definition TBuffer3D.h:113
Box class.
Definition TGeoBBox.h:17
void FillBuffer3D(TBuffer3D &buffer, Int_t reqSections, Bool_t localFrame) const override
Fills the supplied buffer, with sections in desired frame See TBuffer3D.h for explanation of sections...
Double_t fDX
Definition TGeoBBox.h:20
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
Compute distance from outside point to surface of the box.
Definition TGeoBBox.cxx:432
Double_t fOrigin[3]
Definition TGeoBBox.h:23
void InspectShape() const override
Prints shape parameters.
Definition TGeoBBox.cxx:810
Double_t fDY
Definition TGeoBBox.h:21
Double_t fDZ
Definition TGeoBBox.h:22
Int_t GetNsegments() const
Get number of segments approximating circles.
Geometrical transformation package.
Definition TGeoMatrix.h:38
Base abstract class for all shapes.
Definition TGeoShape.h:25
static Double_t Big()
Definition TGeoShape.h:94
Int_t GetBasicColor() const
Get the basic color (0-7).
void TransformPoints(Double_t *points, UInt_t NbPoints) const
Tranform a set of points (LocalToMaster)
void SetShapeBit(UInt_t f, Bool_t set)
Equivalent of TObject::SetBit.
static Double_t SafetyPhi(const Double_t *point, Bool_t in, Double_t phi1, Double_t phi2)
Static method to compute safety w.r.t a phi corner defined by cosines/sines of the angles phi1,...
static Bool_t IsSameWithinTolerance(Double_t a, Double_t b)
Check if two numbers differ with less than a tolerance.
const char * GetPointerName() const
Provide a pointer name containing uid.
Int_t ShapeDistancetoPrimitive(Int_t numpoints, Int_t px, Int_t py) const
Returns distance to shape primitive mesh.
static void NormalPhi(const Double_t *point, const Double_t *dir, Double_t *norm, Double_t c1, Double_t s1, Double_t c2, Double_t s2)
Static method to compute normal to phi planes.
const char * GetName() const override
Get the shape name.
@ kGeoSavePrimitive
Definition TGeoShape.h:64
@ kGeoRunTimeShape
Definition TGeoShape.h:40
static Double_t Tolerance()
Definition TGeoShape.h:97
Bool_t TestShapeBit(UInt_t f) const
Definition TGeoShape.h:175
Int_t SolveQuartic(Double_t a, Double_t b, Double_t c, Double_t d, Double_t *x) const
Find real solutions of the quartic equation : x^4 + a*x^3 + b*x^2 + c*x + d = 0 Input: a,...
void Sizeof3D() const override
fill size of this 3-D object
Double_t ToBoundary(const Double_t *pt, const Double_t *dir, Double_t r, Bool_t in) const
Returns distance to the surface or the torus (fR,r) from a point, along a direction.
Double_t GetAxisRange(Int_t iaxis, Double_t &xlo, Double_t &xhi) const override
Get range of shape for a given axis.
Double_t DDDaxis(const Double_t *pt, const Double_t *dir, Double_t t) const
Second derivative of distance to torus axis w.r.t t.
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
Compute distance from inside point to surface of the torus.
const char * GetAxisName(Int_t iaxis) const override
Returns name of axis IAXIS.
void SavePrimitive(std::ostream &out, Option_t *option="") override
Save a primitive as a C++ statement(s) on output stream "out".
void DistFromInside_v(const Double_t *points, const Double_t *dirs, Double_t *dists, Int_t vecsize, Double_t *step) const override
Compute distance from array of input points having directions specified by dirs. Store output in dist...
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
Compute distance from outside point to surface of the torus.
Double_t Daxis(const Double_t *pt, const Double_t *dir, Double_t t) const
Computes distance to axis of the torus from point pt + t*dir;.
void SetTorusDimensions(Double_t r, Double_t rmin, Double_t rmax, Double_t phi1, Double_t dphi)
Set torus dimensions.
Int_t DistancetoPrimitive(Int_t px, Int_t py) override
Compute closest distance from point px,py to each vertex.
Double_t GetRmin() const
Definition TGeoTorus.h:74
TBuffer3D * MakeBuffer3D() const override
Creates a TBuffer3D describing this shape.
void SetPoints(Double_t *points) const override
Create torus mesh points.
Int_t GetNmeshVertices() const override
Return number of vertices of the mesh representation.
Double_t fR
Definition TGeoTorus.h:20
void SetSegsAndPols(TBuffer3D &buff) const override
Fill TBuffer3D structure for segments and polygons.
void Safety_v(const Double_t *points, const Bool_t *inside, Double_t *safe, Int_t vecsize) const override
Compute safe distance from each of the points in the input array.
TGeoShape * GetMakeRuntimeShape(TGeoShape *mother, TGeoMatrix *mat) const override
Create a shape fitting the mother.
void GetBoundingCylinder(Double_t *param) const override
Fill vector param[4] with the bounding cylinder parameters.
Double_t fRmax
Definition TGeoTorus.h:22
Double_t fDphi
Definition TGeoTorus.h:24
void GetMeshNumbers(Int_t &nvert, Int_t &nsegs, Int_t &npols) const override
Returns numbers of vertices, segments and polygons composing the shape mesh.
Double_t DDaxis(const Double_t *pt, const Double_t *dir, Double_t t) const
Computes derivative w.r.t. t of the distance to axis of the torus from point pt + t*dir;.
Bool_t Contains(const Double_t *point) const override
Test if point is inside the torus.
void Contains_v(const Double_t *points, Bool_t *inside, Int_t vecsize) const override
Check the inside status for each of the points in the array.
void InspectShape() const override
print shape parameters
Double_t fPhi1
Definition TGeoTorus.h:23
TGeoVolume * Divide(TGeoVolume *voldiv, const char *divname, Int_t iaxis, Int_t ndiv, Double_t start, Double_t step) override
Divide this torus shape belonging to volume "voldiv" into ndiv volumes called divname,...
TGeoTorus()
Default constructor.
Definition TGeoTorus.cxx:65
const TBuffer3D & GetBuffer3D(Int_t reqSections, Bool_t localFrame) const override
Fills a static 3D buffer and returns a reference.
void ComputeNormal(const Double_t *point, const Double_t *dir, Double_t *norm) const override
Compute normal to closest surface from POINT.
void SetDimensions(Double_t *param) override
Set torus dimensions starting from a list.
void ComputeBBox() override
Compute bounding box of the torus.
void ComputeNormal_v(const Double_t *points, const Double_t *dirs, Double_t *norms, Int_t vecsize) override
Compute the normal for an array o points so that norm.dot.dir is positive Input: Arrays of point coor...
Int_t SolveCubic(Double_t a, Double_t b, Double_t c, Double_t *x) const
Find real solutions of the cubic equation : x^3 + a*x^2 + b*x + c = 0 Input: a,b,c Output: x[3] real ...
Double_t GetDphi() const
Definition TGeoTorus.h:77
void DistFromOutside_v(const Double_t *points, const Double_t *dirs, Double_t *dists, Int_t vecsize, Double_t *step) const override
Compute distance from array of input points having directions specified by dirs. Store output in dist...
Double_t fRmin
Definition TGeoTorus.h:21
Double_t Safety(const Double_t *point, Bool_t in=kTRUE) const override
computes the closest distance from given point to this shape, according to option.
Double_t Capacity() const override
Computes capacity of the shape in [length^3].
static Double_t DistFromInsideS(const Double_t *point, const Double_t *dir, Double_t rmin, Double_t rmax, Double_t dz, Double_t c1, Double_t s1, Double_t c2, Double_t s2, Double_t cm, Double_t sm, Double_t cdfi)
Compute distance from inside point to surface of the tube segment (static) Boundary safe algorithm.
static Double_t DistFromOutsideS(const Double_t *point, const Double_t *dir, Double_t rmin, Double_t rmax, Double_t dz, Double_t c1, Double_t s1, Double_t c2, Double_t s2, Double_t cm, Double_t sm, Double_t cdfi)
Static method to compute distance to arbitrary tube segment from outside point Boundary safe algorith...
static Double_t DistFromOutsideS(const Double_t *point, const Double_t *dir, Double_t rmin, Double_t rmax, Double_t dz)
Static method to compute distance from outside point to a tube with given parameters Boundary safe al...
Definition TGeoTube.cxx:374
static Double_t DistFromInsideS(const Double_t *point, const Double_t *dir, Double_t rmin, Double_t rmax, Double_t dz)
Compute distance from inside point to surface of the tube (static) Boundary safe algorithm.
Definition TGeoTube.cxx:304
TGeoVolume, TGeoVolumeMulti, TGeoVolumeAssembly are the volume classes.
Definition TGeoVolume.h:43
R__ALWAYS_INLINE Bool_t TestBit(UInt_t f) const
Definition TObject.h:202
virtual const char * ClassName() const
Returns name of class to which the object belongs.
Definition TObject.cxx:226
void SetBit(UInt_t f, Bool_t set)
Set or unset the user status bits as specified in f.
Definition TObject.cxx:864
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition TObject.cxx:1071
TPaveText * pt
Double_t y[n]
Definition legend1.C:17
return c1
Definition legend1.C:41
Double_t x[n]
Definition legend1.C:17
const Int_t n
Definition legend1.C:16
return c2
Definition legend2.C:14
Long64_t LocMin(Long64_t n, const T *a)
Returns index of array with the minimum element.
Definition TMath.h:993
Short_t Max(Short_t a, Short_t b)
Returns the largest of a and b.
Definition TMathBase.h:251
T1 Sign(T1 a, T2 b)
Returns a value with the magnitude of a and the sign of b.
Definition TMathBase.h:176
Double_t ATan2(Double_t y, Double_t x)
Returns the principal value of the arc tangent of y/x, expressed in radians.
Definition TMath.h:657
Long64_t LocMax(Long64_t n, const T *a)
Returns index of array with the maximum element.
Definition TMath.h:1099
constexpr Double_t DegToRad()
Conversion from degree to radian: .
Definition TMath.h:82
Double_t Sqrt(Double_t x)
Returns the square root of x.
Definition TMath.h:673
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Returns x raised to the power y.
Definition TMath.h:732
Short_t Min(Short_t a, Short_t b)
Returns the smallest of a and b.
Definition TMathBase.h:199
Double_t Cos(Double_t)
Returns the cosine of an angle of x radians.
Definition TMath.h:605
constexpr Double_t Pi()
Definition TMath.h:40
Double_t Sin(Double_t)
Returns the sine of an angle of x radians.
Definition TMath.h:599
void Sort(Index n, const Element *a, Index *index, Bool_t down=kTRUE)
Sort the n elements of the array a of generic templated type Element.
Definition TMathBase.h:432
constexpr Double_t RadToDeg()
Conversion from radian to degree: .
Definition TMath.h:75
Short_t Abs(Short_t d)
Returns the absolute value of parameter Short_t d.
Definition TMathBase.h:124