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