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