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