Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
TGeoSphere.cxx
Go to the documentation of this file.
1// @(#)root/geom:$Id$
2// Author: Andrei Gheata 31/01/02
3// TGeoSphere::Contains() DistFromOutside/Out() implemented by Mihaela Gheata
4
5/*************************************************************************
6 * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers. *
7 * All rights reserved. *
8 * *
9 * For the licensing terms see $ROOTSYS/LICENSE. *
10 * For the list of contributors see $ROOTSYS/README/CREDITS. *
11 *************************************************************************/
12
13/** \class TGeoSphere
14\ingroup Shapes_classes
15
16TGeoSphere are not just balls having internal and external
17radii, but sectors of a sphere having defined theta and phi ranges. The
18TGeoSphere class has the following constructor.
19
20~~~{.cpp}
21TGeoSphere(Double_t rmin,Double_t rmax,Double_t theta1,
22Double_t theta2,Double_t phi1, Double_t phi2);
23~~~
24
25Begin_Macro
26{
27 TCanvas *c = new TCanvas("c", "c",0,0,600,600);
28 new TGeoManager("sphere", "poza7");
29 TGeoMaterial *mat = new TGeoMaterial("Al", 26.98,13,2.7);
30 TGeoMedium *med = new TGeoMedium("MED",1,mat);
31 TGeoVolume *top = gGeoManager->MakeBox("TOP",med,100,100,100);
32 gGeoManager->SetTopVolume(top);
33 TGeoVolume *vol = gGeoManager->MakeSphere("SPHERE",med, 30,40,60,120,30,240);
34 vol->SetLineWidth(2);
35 top->AddNode(vol,1);
36 gGeoManager->CloseGeometry();
37 gGeoManager->SetNsegments(30);
38 top->Draw();
39 TView *view = gPad->GetView();
40 if (view) view->ShowAxis();
41}
42End_Macro
43
44 - `rmin: ` internal radius of the spherical sector
45 - `rmax:` external radius
46 - `theta1:` starting theta value [0, 180) in degrees
47 - `theta2:` ending theta value (0, 180] in degrees (`theta1<theta2`)
48
49*/
50
51#include <iostream>
52
53#include "TGeoCone.h"
54#include "TGeoManager.h"
55#include "TGeoVolume.h"
56#include "TVirtualGeoPainter.h"
57#include "TGeoSphere.h"
58#include "TBuffer3D.h"
59#include "TBuffer3DTypes.h"
60#include "TMath.h"
61
62
63////////////////////////////////////////////////////////////////////////////////
64/// Default constructor
65
67{
69 fNz = 0;
70 fNseg = 0;
71 fRmin = 0.0;
72 fRmax = 0.0;
73 fTheta1 = 0.0;
74 fTheta2 = 180.0;
75 fPhi1 = 0.0;
76 fPhi2 = 360.0;
77}
78
79////////////////////////////////////////////////////////////////////////////////
80/// Default constructor specifying minimum and maximum radius
81
90
91////////////////////////////////////////////////////////////////////////////////
92/// Default constructor specifying minimum and maximum radius
93
103
104////////////////////////////////////////////////////////////////////////////////
105/// Default constructor specifying minimum and maximum radius
106/// param[0] = Rmin
107/// param[1] = Rmax
108/// param[2] = theta1
109/// param[3] = theta2
110/// param[4] = phi1
111/// param[5] = phi2
112
120
121////////////////////////////////////////////////////////////////////////////////
122/// destructor
123
125
126////////////////////////////////////////////////////////////////////////////////
127/// Computes capacity of the shape in [length^3]
128
139
140////////////////////////////////////////////////////////////////////////////////
141/// compute bounding box of the sphere
142
144{
148 memset(fOrigin, 0, 3 * sizeof(Double_t));
149 return;
150 }
151 }
159 if (((fTheta1 <= 90) && (fTheta2 >= 90)) || ((fTheta2 <= 90) && (fTheta1 >= 90))) {
160 r1max = fRmax;
161 r2max = fRmin;
162 }
165
166 Double_t xc[4];
167 Double_t yc[4];
176
177 Double_t xmin = xc[TMath::LocMin(4, &xc[0])];
178 Double_t xmax = xc[TMath::LocMax(4, &xc[0])];
179 Double_t ymin = yc[TMath::LocMin(4, &yc[0])];
180 Double_t ymax = yc[TMath::LocMax(4, &yc[0])];
182 if (dp < 0)
183 dp += 360;
184 Double_t ddp = -fPhi1;
185 if (ddp < 0)
186 ddp += 360;
187 if (ddp > 360)
188 ddp -= 360;
189 if (ddp <= dp)
190 xmax = rmax;
191 ddp = 90 - fPhi1;
192 if (ddp < 0)
193 ddp += 360;
194 if (ddp > 360)
195 ddp -= 360;
196 if (ddp <= dp)
197 ymax = rmax;
198 ddp = 180 - fPhi1;
199 if (ddp < 0)
200 ddp += 360;
201 if (ddp > 360)
202 ddp -= 360;
203 if (ddp <= dp)
204 xmin = -rmax;
205 ddp = 270 - fPhi1;
206 if (ddp < 0)
207 ddp += 360;
208 if (ddp > 360)
209 ddp -= 360;
210 if (ddp <= dp)
211 ymin = -rmax;
216 Double_t zmin = xc[TMath::LocMin(4, &xc[0])];
217 Double_t zmax = xc[TMath::LocMax(4, &xc[0])];
218
219 fOrigin[0] = (xmax + xmin) / 2;
220 fOrigin[1] = (ymax + ymin) / 2;
221 fOrigin[2] = (zmax + zmin) / 2;
222 ;
223 fDX = (xmax - xmin) / 2;
224 fDY = (ymax - ymin) / 2;
225 fDZ = (zmax - zmin) / 2;
226}
227
228////////////////////////////////////////////////////////////////////////////////
229/// Compute normal to closest surface from POINT.
230
231void TGeoSphere::ComputeNormal(const Double_t *point, const Double_t *dir, Double_t *norm) const
232{
233 Double_t rxy2 = point[0] * point[0] + point[1] * point[1];
234 Double_t r2 = rxy2 + point[2] * point[2];
237 if (r <= 1E-20)
238 rzero = kTRUE;
239 // localize theta
240 Double_t phi = 0;
241 Double_t th = 0.;
242 if (!rzero)
243 th = TMath::ACos(point[2] / r);
244
245 // localize phi
246 phi = TMath::ATan2(point[1], point[0]);
247
248 Double_t saf[4];
251 : TMath::Abs(r - fRmin);
252 saf[1] = TMath::Abs(fRmax - r);
253 saf[2] = saf[3] = TGeoShape::Big();
255 if (fTheta1 > 0) {
257 }
258 if (fTheta2 < 180) {
260 }
261 }
262 Int_t i = TMath::LocMin(4, saf);
268 if (TGeoShape::IsCloseToPhi(saf[i], point, c1, s1, c2, s2)) {
269 TGeoShape::NormalPhi(point, dir, norm, c1, s1, c2, s2);
270 return;
271 }
272 }
273 if (i > 1) {
274 if (i == 2)
275 th = (fTheta1 < 90) ? (fTheta1 + 90) : (fTheta1 - 90);
276 else
277 th = (fTheta2 < 90) ? (fTheta2 + 90) : (fTheta2 - 90);
278 th *= TMath::DegToRad();
279 }
280
281 norm[0] = TMath::Sin(th) * TMath::Cos(phi);
282 norm[1] = TMath::Sin(th) * TMath::Sin(phi);
283 norm[2] = TMath::Cos(th);
284 if (norm[0] * dir[0] + norm[1] * dir[1] + norm[2] * dir[2] < 0) {
285 norm[0] = -norm[0];
286 norm[1] = -norm[1];
287 norm[2] = -norm[2];
288 }
289}
290
291////////////////////////////////////////////////////////////////////////////////
292/// Check if a point in local sphere coordinates is close to a boundary within
293/// shape tolerance. Return values:
294/// - 0 - not close to boundary
295/// - 1 - close to Rmin boundary
296/// - 2 - close to Rmax boundary
297/// - 3,4 - close to phi1/phi2 boundary
298/// - 5,6 - close to theta1/theta2 boundary
299
301{
302 Int_t icode = 0;
304 Double_t r2 = point[0] * point[0] + point[1] * point[1] + point[2] * point[2];
306 // Test if point is on fRmax boundary
307 if (TMath::Abs(drsqout) < 2. * fRmax * tol)
308 return 2;
310 // Test if point is on fRmin boundary
311 if (TestShapeBit(kGeoRSeg)) {
312 drsqin -= fRmin * fRmin;
313 if (TMath::Abs(drsqin) < 2. * fRmin * tol)
314 return 1;
315 }
317 Double_t phi = TMath::ATan2(point[1], point[0]);
318 if (phi < 0)
319 phi += 2 * TMath::Pi();
322 Double_t ddp = phi - phi1;
323 if (r2 * ddp * ddp < tol * tol)
324 return 3;
325 ddp = phi - phi2;
326 if (r2 * ddp * ddp < tol * tol)
327 return 4;
328 }
331 Double_t theta = TMath::ACos(point[2] / r2);
335 if (fTheta1 > 0) {
336 ddt = TMath::Abs(theta - theta1);
337 if (r * ddt < tol)
338 return 5;
339 }
340 if (fTheta2 < 180) {
341 ddt = TMath::Abs(theta - theta2);
342 if (r * ddt < tol)
343 return 6;
344 }
345 }
346 return icode;
347}
348
349////////////////////////////////////////////////////////////////////////////////
350/// Check if a point is inside radius/theta/phi ranges for the spherical sector.
351
353{
354 Double_t r2 = point[0] * point[0] + point[1] * point[1] + point[2] * point[2];
355 if (checkR) {
356 if (TestShapeBit(kGeoRSeg) && (r2 < fRmin * fRmin))
357 return kFALSE;
358 if (r2 > fRmax * fRmax)
359 return kFALSE;
360 }
361 if (r2 < 1E-20)
362 return kTRUE;
364 Double_t phi = TMath::ATan2(point[1], point[0]) * TMath::RadToDeg();
365 while (phi < fPhi1)
366 phi += 360.;
368 Double_t ddp = phi - fPhi1;
369 if (ddp > dphi)
370 return kFALSE;
371 }
373 r2 = TMath::Sqrt(r2);
374 // check theta range
375 Double_t theta = TMath::ACos(point[2] / r2) * TMath::RadToDeg();
376 if ((theta < fTheta1) || (theta > fTheta2))
377 return kFALSE;
378 }
379 return kTRUE;
380}
381
382////////////////////////////////////////////////////////////////////////////////
383/// test if point is inside this sphere
384/// check Rmin<=R<=Rmax
385
387{
388 Double_t r2 = point[0] * point[0] + point[1] * point[1] + point[2] * point[2];
389 if (TestShapeBit(kGeoRSeg) && (r2 < fRmin * fRmin))
390 return kFALSE;
391 if (r2 > fRmax * fRmax)
392 return kFALSE;
393 if (r2 < 1E-20)
394 return kTRUE;
395 // check phi range
397 Double_t phi = TMath::ATan2(point[1], point[0]) * TMath::RadToDeg();
398 if (phi < 0)
399 phi += 360.;
401 if (dphi < 0)
402 dphi += 360.;
403 Double_t ddp = phi - fPhi1;
404 if (ddp < 0)
405 ddp += 360.;
406 if (ddp > dphi)
407 return kFALSE;
408 }
410 r2 = TMath::Sqrt(r2);
411 // check theta range
412 Double_t theta = TMath::ACos(point[2] / r2) * TMath::RadToDeg();
413 if ((theta < fTheta1) || (theta > fTheta2))
414 return kFALSE;
415 }
416 return kTRUE;
417}
418
419////////////////////////////////////////////////////////////////////////////////
420/// compute closest distance from point px,py to each corner
421
423{
424 Int_t n = fNseg + 1;
425 Int_t nz = fNz + 1;
426 const Int_t numPoints = 2 * n * nz;
427 return ShapeDistancetoPrimitive(numPoints, px, py);
428}
429
430////////////////////////////////////////////////////////////////////////////////
431/// compute distance from outside point to surface of the sphere
432/// Check if the bounding box is crossed within the requested distance
433
436{
437 Double_t sdist = TGeoBBox::DistFromOutside(point, dir, fDX, fDY, fDZ, fOrigin, step);
438 if (sdist >= step)
439 return TGeoShape::Big();
440 if (iact < 3 && safe) {
441 // Reuse the outside lower bound, including the angular-cut constraints.
442 *safe = Safety(point, kFALSE);
443 if (iact == 0)
444 return TGeoShape::Big();
445 if (iact == 1 && step < *safe)
446 return TGeoShape::Big();
447 }
448 // compute distance to shape
449 Double_t r1, r2, z1, z2, dz, si, ci;
450 Double_t rxy2 = point[0] * point[0] + point[1] * point[1];
451 r2 = rxy2 + point[2] * point[2];
453 // first check if any crossing at all
455 Double_t rdotn = point[0] * dir[0] + point[1] * dir[1] + point[2] * dir[2];
457 if (r > fRmax) {
458 Double_t b = rdotn;
459 Double_t c = r2 - fRmax * fRmax;
460 Double_t d = b * b - c;
461 if (d < 0)
462 return TGeoShape::Big();
463 }
464 if (fullsph) {
467 if (r <= fRmax + TGeoShape::Tolerance())
468 inrmax = kTRUE;
469 if (r >= fRmin - TGeoShape::Tolerance())
470 inrmin = kTRUE;
471 if (inrmax && inrmin) {
472 if ((fRmax - r) < (r - fRmin)) {
473 // close to Rmax
474 if (rdotn >= 0)
475 return TGeoShape::Big();
476 return 0.0; // already in
477 }
478 // close to Rmin
480 return 0.0;
481 // check second crossing of Rmin
482 return DistToSphere(point, dir, fRmin, kFALSE, kFALSE);
483 }
484 }
485
486 // do rmin, rmax, checking phi and theta ranges
487 if (r < fRmin) {
488 // check first cross of rmin
489 snxt = DistToSphere(point, dir, fRmin, kTRUE);
490 if (snxt < 1E20)
491 return snxt;
492 } else {
493 if (r > fRmax) {
494 // point outside rmax, check first cross of rmax
495 snxt = DistToSphere(point, dir, fRmax, kTRUE);
496 if (snxt < 1E20)
497 return snxt;
498 // now check second crossing of rmin
499 if (fRmin > 0)
500 snxt = DistToSphere(point, dir, fRmin, kTRUE, kFALSE);
501 } else {
502 // point between rmin and rmax, check second cross of rmin
503 if (fRmin > 0)
504 snxt = DistToSphere(point, dir, fRmin, kTRUE, kFALSE);
505 }
506 }
507 // check theta conical surfaces
508 Double_t ptnew[3];
509 Double_t b, delta, xnew, ynew, znew, phi0, ddp;
512
514 if (fTheta1 > 0) {
516 // surface is a plane
517 if (point[2] * dir[2] < 0) {
518 snxt = -point[2] / dir[2];
519 ptnew[0] = point[0] + snxt * dir[0];
520 ptnew[1] = point[1] + snxt * dir[1];
521 ptnew[2] = 0;
522 // check range
524 return TMath::Min(snxt, snext);
525 }
526 } else {
529 if (ci > 0) {
530 r1 = fRmin * si;
531 z1 = fRmin * ci;
532 r2 = fRmax * si;
533 z2 = fRmax * ci;
534 } else {
535 r1 = fRmax * si;
536 z1 = fRmax * ci;
537 r2 = fRmin * si;
538 z2 = fRmin * ci;
539 }
540 dz = 0.5 * (z2 - z1);
541 ptnew[0] = point[0];
542 ptnew[1] = point[1];
543 ptnew[2] = point[2] - 0.5 * (z1 + z2);
544 Double_t zinv = 1. / dz;
545 Double_t rin = 0.5 * (r1 + r2 + (r2 - r1) * ptnew[2] * zinv);
546 // Protection in case point is outside
547 Double_t sigz = TMath::Sign(1., point[2]);
548 TGeoCone::DistToCone(ptnew, dir, dz, r1, r2, b, delta);
550 if (delta < 0)
551 skip = kTRUE;
552 if (!skip) {
553 snxt = -b - delta;
554 if (sigz * ci > 0 && sigz * rxy2 > sigz * rin * (rin - sigz * TGeoShape::Tolerance())) {
556 ptnew[0] * dir[0] + ptnew[1] * dir[1] + 0.5 * (r1 - r2) * dir[2] * zinv * TMath::Sqrt(rxy2);
557 if (sigz * ddotn >= 0 || -b + delta < 1.E-9)
558 skip = kTRUE;
560 }
561 if (snxt < 1E10) {
562 znew = ptnew[2] + snxt * dir[2];
563 if (snxt > 0 && TMath::Abs(znew) < dz) {
565 st1 = snxt;
566 else {
567 xnew = ptnew[0] + snxt * dir[0];
568 ynew = ptnew[1] + snxt * dir[1];
570 ddp = phi0 - fPhi1;
571 while (ddp < 0)
572 ddp += 360.;
573 if (ddp <= fPhi2 - fPhi1)
574 st1 = snxt;
575 }
576 }
577 }
578 if (!skip && st1 > 1E10) {
579 snxt = -b + delta;
580 znew = ptnew[2] + snxt * dir[2];
581 if (snxt > 0 && TMath::Abs(znew) < dz) {
583 st1 = snxt;
584 else {
585 xnew = ptnew[0] + snxt * dir[0];
586 ynew = ptnew[1] + snxt * dir[1];
588 ddp = phi0 - fPhi1;
589 while (ddp < 0)
590 ddp += 360.;
591 if (ddp <= fPhi2 - fPhi1)
592 st1 = snxt;
593 }
594 }
595 }
596 }
597 }
598 }
599
600 if (fTheta2 < 180) {
602 // surface is a plane
603 if (point[2] * dir[2] < 0) {
604 snxt = -point[2] / dir[2];
605 ptnew[0] = point[0] + snxt * dir[0];
606 ptnew[1] = point[1] + snxt * dir[1];
607 ptnew[2] = 0;
608 // check range
610 return TMath::Min(snxt, snext);
611 }
612 } else {
615 if (ci > 0) {
616 r1 = fRmin * si;
617 z1 = fRmin * ci;
618 r2 = fRmax * si;
619 z2 = fRmax * ci;
620 } else {
621 r1 = fRmax * si;
622 z1 = fRmax * ci;
623 r2 = fRmin * si;
624 z2 = fRmin * ci;
625 }
626 dz = 0.5 * (z2 - z1);
627 ptnew[0] = point[0];
628 ptnew[1] = point[1];
629 ptnew[2] = point[2] - 0.5 * (z1 + z2);
630 Double_t zinv = 1. / dz;
631 Double_t rin = 0.5 * (r1 + r2 + (r2 - r1) * ptnew[2] * zinv);
632 // Protection in case point is outside
633 Double_t sigz = TMath::Sign(1., point[2]);
634 TGeoCone::DistToCone(ptnew, dir, dz, r1, r2, b, delta);
636 if (delta < 0)
637 skip = kTRUE;
638 if (!skip) {
639 snxt = -b - delta;
640 if (sigz * ci > 0 && sigz * rxy2 < sigz * rin * (rin + sigz * TGeoShape::Tolerance())) {
642 ptnew[0] * dir[0] + ptnew[1] * dir[1] + 0.5 * (r1 - r2) * dir[2] * zinv * TMath::Sqrt(rxy2);
643 if (sigz * ddotn <= 0 || -b + delta < 1.E-9)
644 skip = kTRUE;
646 }
647 if (snxt < 1E10) {
648 znew = ptnew[2] + snxt * dir[2];
649 if (snxt > 0 && TMath::Abs(znew) < dz) {
651 st2 = snxt;
652 else {
653 xnew = ptnew[0] + snxt * dir[0];
654 ynew = ptnew[1] + snxt * dir[1];
656 ddp = phi0 - fPhi1;
657 while (ddp < 0)
658 ddp += 360.;
659 if (ddp <= fPhi2 - fPhi1)
660 st2 = snxt;
661 }
662 }
663 }
664 if (!skip && st2 > 1E10) {
665 snxt = -b + delta;
666 znew = ptnew[2] + snxt * dir[2];
667 if (snxt > 0 && TMath::Abs(znew) < dz) {
669 st2 = snxt;
670 else {
671 xnew = ptnew[0] + snxt * dir[0];
672 ynew = ptnew[1] + snxt * dir[1];
674 ddp = phi0 - fPhi1;
675 while (ddp < 0)
676 ddp += 360.;
677 if (ddp <= fPhi2 - fPhi1)
678 st2 = snxt;
679 }
680 }
681 }
682 }
683 }
684 }
685 }
688 // if (snxt<1E20) return snxt;
694 Double_t phim = 0.5 * (fPhi1 + fPhi2);
697 Double_t s = 0;
699 safety = point[0] * s1 - point[1] * c1;
700 if (safety > 0) {
701 un = dir[0] * s1 - dir[1] * c1;
702 if (un < 0) {
703 s = -safety / un;
704 ptnew[0] = point[0] + s * dir[0];
705 ptnew[1] = point[1] + s * dir[1];
706 ptnew[2] = point[2] + s * dir[2];
707 if ((ptnew[1] * cm - ptnew[0] * sm) <= 0) {
708 Double_t sfi1 = s;
709 if (IsPointInside(&ptnew[0], kTRUE, kTRUE, kFALSE) && sfi1 < snxt)
710 return sfi1;
711 }
712 }
713 }
714 safety = -point[0] * s2 + point[1] * c2;
715 if (safety > 0) {
716 un = -dir[0] * s2 + dir[1] * c2;
717 if (un < 0) {
718 s = -safety / un;
719 ptnew[0] = point[0] + s * dir[0];
720 ptnew[1] = point[1] + s * dir[1];
721 ptnew[2] = point[2] + s * dir[2];
722 if ((ptnew[1] * cm - ptnew[0] * sm) >= 0) {
723 Double_t sfi2 = s;
724 if (IsPointInside(&ptnew[0], kTRUE, kTRUE, kFALSE) && sfi2 < snxt)
725 return sfi2;
726 }
727 }
728 }
729 }
730 return snxt;
731}
732
733////////////////////////////////////////////////////////////////////////////////
734/// compute distance from inside point to surface of the sphere
735
738{
739 Double_t saf[6];
740 Double_t rxy2 = point[0] * point[0] + point[1] * point[1];
742 Double_t rad2 = rxy2 + point[2] * point[2];
745 if (r <= 1E-20)
746 rzero = kTRUE;
747 // localize theta
748 Double_t phi = 0;
749 ;
750 Double_t th = 0.;
751 if (TestShapeBit(kGeoThetaSeg) && (!rzero)) {
752 th = TMath::ACos(point[2] / r) * TMath::RadToDeg();
753 }
754 // localize phi
756 phi = TMath::ATan2(point[1], point[0]) * TMath::RadToDeg();
757 if (phi < 0)
758 phi += 360.;
759 }
760 if (iact < 3 && safe) {
762 saf[1] = fRmax - r;
763 saf[2] = saf[3] = saf[4] = saf[5] = TGeoShape::Big();
765 if (fTheta1 > 0) {
766 saf[2] = r * TMath::Sin((th - fTheta1) * TMath::DegToRad());
767 }
768 if (fTheta2 < 180) {
769 saf[3] = r * TMath::Sin((fTheta2 - th) * TMath::DegToRad());
770 }
771 }
773 Double_t dph1 = phi - fPhi1;
774 if (dph1 < 0)
775 dph1 += 360.;
776 if (dph1 <= 90.)
778 Double_t dph2 = fPhi2 - phi;
779 if (dph2 < 0)
780 dph2 += 360.;
781 if (dph2 <= 90.)
783 }
784 *safe = saf[TMath::LocMin(6, &saf[0])];
785 if (iact == 0)
786 return TGeoShape::Big();
787 if (iact == 1 && step < *safe)
788 return TGeoShape::Big();
789 }
790 // compute distance to shape
791 if (rzero) {
792 // gGeoManager->SetNormalChecked(1.);
793 return fRmax;
794 }
795 // first do rmin, rmax
796 Double_t b, delta, xnew, ynew, znew, phi0, ddp;
797 Double_t rdotn = point[0] * dir[0] + point[1] * dir[1] + point[2] * dir[2];
799 // Inner sphere
800 if (fRmin > 0) {
801 // Protection in case point is actually outside the sphere
802 if (r <= fRmin + TGeoShape::Tolerance()) {
803 if (rdotn < 0)
804 return 0.0;
805 } else {
806 if (rdotn < 0)
807 sn1 = DistToSphere(point, dir, fRmin, kFALSE);
808 }
809 }
810 // Outer sphere
811 if (r >= fRmax - TGeoShape::Tolerance()) {
812 if (rdotn >= 0)
813 return 0.0;
814 }
815 Double_t sn2 = DistToSphere(point, dir, fRmax, kFALSE, kFALSE);
816 Double_t sr = TMath::Min(sn1, sn2);
817 // check theta conical surfaces
818 sn1 = sn2 = TGeoShape::Big();
821 // surface is a plane
822 if (point[2] * dir[2] < 0)
823 sn1 = -point[2] / dir[2];
824 } else {
825 if (fTheta1 > 0) {
826 Double_t r1, r2, z1, z2, dz, ptnew[3];
829 if (ci > 0) {
830 r1 = fRmin * si;
831 z1 = fRmin * ci;
832 r2 = fRmax * si;
833 z2 = fRmax * ci;
834 } else {
835 r1 = fRmax * si;
836 z1 = fRmax * ci;
837 r2 = fRmin * si;
838 z2 = fRmin * ci;
839 }
840 dz = 0.5 * (z2 - z1);
841 ptnew[0] = point[0];
842 ptnew[1] = point[1];
843 ptnew[2] = point[2] - 0.5 * (z1 + z2);
844 Double_t zinv = 1. / dz;
845 Double_t rin = 0.5 * (r1 + r2 + (r2 - r1) * ptnew[2] * zinv);
846 // Protection in case point is outside
847 Double_t sigz = TMath::Sign(1., point[2]);
848 if (sigz * ci > 0 && sigz * rxy2 < sigz * rin * (rin + sigz * TGeoShape::Tolerance())) {
850 ptnew[0] * dir[0] + ptnew[1] * dir[1] + 0.5 * (r1 - r2) * dir[2] * zinv * TMath::Sqrt(rxy2);
851 if (sigz * ddotn <= 0)
852 return 0.0;
853 } else {
854 TGeoCone::DistToCone(ptnew, dir, dz, r1, r2, b, delta);
855 if (delta > 0) {
856 Double_t snxt = -b - delta;
857 znew = ptnew[2] + snxt * dir[2];
858 if (snxt > 0 && TMath::Abs(znew) < dz) {
860 sn1 = snxt;
861 else {
862 xnew = ptnew[0] + snxt * dir[0];
863 ynew = ptnew[1] + snxt * dir[1];
865 ddp = phi0 - fPhi1;
866 while (ddp < 0)
867 ddp += 360.;
868 if (ddp <= fPhi2 - fPhi1)
869 sn1 = snxt;
870 }
871 }
872 if (sn1 > 1E10) {
873 snxt = -b + delta;
874 znew = ptnew[2] + snxt * dir[2];
875 if (snxt > 0 && TMath::Abs(znew) < dz) {
877 sn1 = snxt;
878 else {
879 xnew = ptnew[0] + snxt * dir[0];
880 ynew = ptnew[1] + snxt * dir[1];
882 ddp = phi0 - fPhi1;
883 while (ddp < 0)
884 ddp += 360.;
885 if (ddp <= fPhi2 - fPhi1)
886 sn1 = snxt;
887 }
888 }
889 }
890 }
891 }
892 }
893 }
895 // surface is a plane
896 if (point[2] * dir[2] < 0)
897 sn1 = -point[2] / dir[2];
898 } else {
899 if (fTheta2 < 180) {
900 Double_t r1, r2, z1, z2, dz, ptnew[3];
903 if (ci > 0) {
904 r1 = fRmin * si;
905 z1 = fRmin * ci;
906 r2 = fRmax * si;
907 z2 = fRmax * ci;
908 } else {
909 r1 = fRmax * si;
910 z1 = fRmax * ci;
911 r2 = fRmin * si;
912 z2 = fRmin * ci;
913 }
914 dz = 0.5 * (z2 - z1);
915 ptnew[0] = point[0];
916 ptnew[1] = point[1];
917 ptnew[2] = point[2] - 0.5 * (z1 + z2);
918 Double_t zinv = 1. / dz;
919 Double_t rin = 0.5 * (r1 + r2 + (r2 - r1) * ptnew[2] * zinv);
920 // Protection in case point is outside
921 Double_t sigz = TMath::Sign(1., point[2]);
922 if (sigz * ci > 0 && sigz * rxy2 > sigz * rin * (rin - sigz * TGeoShape::Tolerance())) {
924 ptnew[0] * dir[0] + ptnew[1] * dir[1] + 0.5 * (r1 - r2) * dir[2] * zinv * TMath::Sqrt(rxy2);
925 if (sigz * ddotn >= 0)
926 return 0.0;
927 } else {
928 TGeoCone::DistToCone(ptnew, dir, dz, r1, r2, b, delta);
929 if (delta > 0) {
930 Double_t snxt = -b - delta;
931 znew = ptnew[2] + snxt * dir[2];
932 if (snxt > 0 && TMath::Abs(znew) < dz) {
934 sn2 = snxt;
935 else {
936 xnew = ptnew[0] + snxt * dir[0];
937 ynew = ptnew[1] + snxt * dir[1];
939 ddp = phi0 - fPhi1;
940 while (ddp < 0)
941 ddp += 360.;
942 if (ddp <= fPhi2 - fPhi1)
943 sn2 = snxt;
944 }
945 }
946 if (sn2 > 1E10) {
947 snxt = -b + delta;
948 znew = ptnew[2] + snxt * dir[2];
949 if (snxt > 0 && TMath::Abs(znew) < dz) {
951 sn2 = snxt;
952 else {
953 xnew = ptnew[0] + snxt * dir[0];
954 ynew = ptnew[1] + snxt * dir[1];
956 ddp = phi0 - fPhi1;
957 while (ddp < 0)
958 ddp += 360.;
959 if (ddp <= fPhi2 - fPhi1)
960 sn2 = snxt;
961 }
962 }
963 }
964 }
965 }
966 }
967 }
968 }
976 Double_t phim = 0.5 * (fPhi1 + fPhi2);
979 sp = TGeoShape::DistToPhiMin(point, dir, s1, c1, s2, c2, sm, cm);
980 }
981 Double_t snxt = TMath::Min(sr, st);
983 return snxt;
984}
985
986////////////////////////////////////////////////////////////////////////////////
987/// compute distance to sphere of radius rsph. Direction has to be a unit vector
988
990 Bool_t firstcross) const
991{
992 if (rsph <= 0)
993 return TGeoShape::Big();
994 Double_t r2 = point[0] * point[0] + point[1] * point[1] + point[2] * point[2];
995 Double_t b = point[0] * dir[0] + point[1] * dir[1] + point[2] * dir[2];
996 Double_t c = r2 - rsph * rsph;
997 Bool_t in = (c <= 0) ? kTRUE : kFALSE;
998 Double_t d;
999
1000 d = b * b - c;
1001 if (d < 0)
1002 return TGeoShape::Big();
1003 Double_t pt[3];
1004 Int_t i;
1005 d = TMath::Sqrt(d);
1006 Double_t s;
1007 if (in) {
1008 s = -b + d;
1009 } else {
1010 s = (firstcross) ? (-b - d) : (-b + d);
1011 }
1012 if (s < 0)
1013 return TGeoShape::Big();
1014 if (!check)
1015 return s;
1016 for (i = 0; i < 3; i++)
1017 pt[i] = point[i] + s * dir[i];
1018 // check theta and phi ranges
1019 if (IsPointInside(&pt[0], kFALSE))
1020 return s;
1021 return TGeoShape::Big();
1022}
1023
1024////////////////////////////////////////////////////////////////////////////////
1025
1026TGeoVolume *
1028{
1029 TGeoShape *shape; //--- shape to be created
1030 TGeoVolume *vol; //--- division volume to be created
1031 TGeoVolumeMulti *vmulti; //--- generic divided volume
1032 TGeoPatternFinder *finder; //--- finder to be attached
1033 TString opt = ""; //--- option to be attached
1034 Int_t id;
1035 Double_t end = start + ndiv * step;
1036 switch (iaxis) {
1037 case 1: //--- R division
1038 finder = new TGeoPatternSphR(voldiv, ndiv, start, end);
1040 voldiv->SetFinder(finder);
1041 finder->SetDivIndex(voldiv->GetNdaughters());
1042 for (id = 0; id < ndiv; id++) {
1043 shape = new TGeoSphere(start + id * step, start + (id + 1) * step, fTheta1, fTheta2, fPhi1, fPhi2);
1044 vol = new TGeoVolume(divname, shape, voldiv->GetMedium());
1045 vmulti->AddVolume(vol);
1046 opt = "R";
1047 voldiv->AddNodeOffset(vol, id, 0, opt.Data());
1048 ((TGeoNodeOffset *)voldiv->GetNodes()->At(voldiv->GetNdaughters() - 1))->SetFinder(finder);
1049 }
1050 return vmulti;
1051 case 2: //--- Phi division
1052 finder = new TGeoPatternSphPhi(voldiv, ndiv, start, end);
1053 voldiv->SetFinder(finder);
1054 finder->SetDivIndex(voldiv->GetNdaughters());
1055 shape = new TGeoSphere(fRmin, fRmax, fTheta1, fTheta2, -step / 2, step / 2);
1056 vol = new TGeoVolume(divname, shape, voldiv->GetMedium());
1058 vmulti->AddVolume(vol);
1059 opt = "Phi";
1060 for (id = 0; id < ndiv; id++) {
1061 voldiv->AddNodeOffset(vol, id, start + id * step + step / 2, opt.Data());
1062 ((TGeoNodeOffset *)voldiv->GetNodes()->At(voldiv->GetNdaughters() - 1))->SetFinder(finder);
1063 }
1064 return vmulti;
1065 case 3: //--- Theta division
1066 finder = new TGeoPatternSphTheta(voldiv, ndiv, start, end);
1068 voldiv->SetFinder(finder);
1069 finder->SetDivIndex(voldiv->GetNdaughters());
1070 for (id = 0; id < ndiv; id++) {
1071 shape = new TGeoSphere(fRmin, fRmax, start + id * step, start + (id + 1) * step, fPhi1, fPhi2);
1072 vol = new TGeoVolume(divname, shape, voldiv->GetMedium());
1073 vmulti->AddVolume(vol);
1074 opt = "Theta";
1075 voldiv->AddNodeOffset(vol, id, 0, opt.Data());
1076 ((TGeoNodeOffset *)voldiv->GetNodes()->At(voldiv->GetNdaughters() - 1))->SetFinder(finder);
1077 }
1078 return vmulti;
1079 default: Error("Divide", "In shape %s wrong axis type for division", GetName()); return nullptr;
1080 }
1081}
1082
1083////////////////////////////////////////////////////////////////////////////////
1084/// Returns name of axis IAXIS.
1085
1087{
1088 switch (iaxis) {
1089 case 1: return "R";
1090 case 2: return "PHI";
1091 case 3: return "THETA";
1092 default: return "UNDEFINED";
1093 }
1094}
1095
1096////////////////////////////////////////////////////////////////////////////////
1097/// Get range of shape for a given axis.
1098
1100{
1101 xlo = 0;
1102 xhi = 0;
1103 Double_t dx = 0;
1104 switch (iaxis) {
1105 case 1:
1106 xlo = fRmin;
1107 xhi = fRmax;
1108 dx = xhi - xlo;
1109 return dx;
1110 case 2:
1111 xlo = fPhi1;
1112 xhi = fPhi2;
1113 dx = xhi - xlo;
1114 return dx;
1115 case 3:
1116 xlo = fTheta1;
1117 xhi = fTheta2;
1118 dx = xhi - xlo;
1119 return dx;
1120 }
1121 return dx;
1122}
1123
1124////////////////////////////////////////////////////////////////////////////////
1125/// Fill vector param[4] with the bounding cylinder parameters. The order
1126/// is the following : Rmin, Rmax, Phi1, Phi2
1127
1129{
1132 if (smin > smax) {
1133 Double_t a = smin;
1134 smin = smax;
1135 smax = a;
1136 }
1137 param[0] = fRmin * smin; // Rmin
1138 param[0] *= param[0];
1139 if (((90. - fTheta1) * (fTheta2 - 90.)) >= 0)
1140 smax = 1.;
1141 param[1] = fRmax * smax; // Rmax
1142 param[1] *= param[1];
1143 param[2] = (fPhi1 < 0) ? (fPhi1 + 360.) : fPhi1; // Phi1
1144 param[3] = fPhi2;
1145 if (TGeoShape::IsSameWithinTolerance(param[3] - param[2], 360)) { // Phi2
1146 param[2] = 0.;
1147 param[3] = 360.;
1148 }
1149 while (param[3] < param[2])
1150 param[3] += 360.;
1151}
1152
1153////////////////////////////////////////////////////////////////////////////////
1154/// print shape parameters
1155
1157{
1158 printf("*** Shape %s: TGeoSphere ***\n", GetName());
1159 printf(" Rmin = %11.5f\n", fRmin);
1160 printf(" Rmax = %11.5f\n", fRmax);
1161 printf(" Th1 = %11.5f\n", fTheta1);
1162 printf(" Th2 = %11.5f\n", fTheta2);
1163 printf(" Ph1 = %11.5f\n", fPhi1);
1164 printf(" Ph2 = %11.5f\n", fPhi2);
1165 printf(" Bounding box:\n");
1167}
1168
1169////////////////////////////////////////////////////////////////////////////////
1170/// Creates a TBuffer3D describing *this* shape.
1171/// Coordinates are in local reference frame.
1172
1174{
1175 Bool_t full = kTRUE;
1177 full = kFALSE;
1178 Int_t ncenter = 1;
1179 if (full || TestShapeBit(kGeoRSeg))
1180 ncenter = 0;
1181 Int_t nup = (fTheta1 > 0) ? 0 : 1;
1182 Int_t ndown = (fTheta2 < 180) ? 0 : 1;
1183 // number of different latitudes, excluding 0 and 180 degrees
1184 Int_t nlat = fNz + 1 - (nup + ndown);
1185 // number of different longitudes
1186 Int_t nlong = fNseg;
1188 nlong++;
1189
1192 nbPnts *= 2;
1193
1194 Int_t nbSegs = nlat * fNseg + (nlat - 1 + nup + ndown) * nlong; // outer sphere
1196 nbSegs *= 2; // inner sphere
1198 nbSegs += 2 * nlat + nup + ndown; // 2 phi planes
1199 nbSegs += nlong * (2 - nup - ndown); // connecting cones
1200
1201 Int_t nbPols = fNz * fNseg; // outer
1203 nbPols *= 2; // inner
1205 nbPols += 2 * fNz; // 2 phi planes
1206 nbPols += (2 - nup - ndown) * fNseg; // connecting
1207
1208 TBuffer3D *buff =
1210
1211 if (buff) {
1212 SetPoints(buff->fPnts);
1214 }
1215
1216 return buff;
1217}
1218
1219////////////////////////////////////////////////////////////////////////////////
1220/// Fill TBuffer3D structure for segments and polygons.
1221
1223{
1224 // Bool_t full = kTRUE;
1225 // if (TestShapeBit(kGeoThetaSeg) || TestShapeBit(kGeoPhiSeg)) full = kFALSE;
1226 // Int_t ncenter = 1;
1227 // if (full || TestShapeBit(kGeoRSeg)) ncenter = 0;
1228 Int_t nup = (fTheta1 > 0) ? 0 : 1;
1229 Int_t ndown = (fTheta2 < 180) ? 0 : 1;
1230 // number of different latitudes, excluding 0 and 180 degrees
1231 Int_t nlat = fNz + 1 - (nup + ndown);
1232 // number of different longitudes
1233 Int_t nlong = fNseg;
1235 nlong++;
1236
1237 // Int_t nbPnts = nlat*nlong+nup+ndown+ncenter;
1238 // if (TestShapeBit(kGeoRSeg)) nbPnts *= 2;
1239
1240 // Int_t nbSegs = nlat*fNseg + (nlat-1+nup+ndown)*nlong; // outer sphere
1241 // if (TestShapeBit(kGeoRSeg)) nbSegs *= 2; // inner sphere
1242 // if (TestShapeBit(kGeoPhiSeg)) nbSegs += 2*nlat+nup+ndown; // 2 phi planes
1243 // nbSegs += nlong * (2-nup - ndown); // connecting cones
1244
1245 // Int_t nbPols = fNz*fNseg; // outer
1246 // if (TestShapeBit(kGeoRSeg)) nbPols *=2; // inner
1247 // if (TestShapeBit(kGeoPhiSeg)) nbPols += 2*fNz; // 2 phi planes
1248 // nbPols += (2-nup-ndown)*fNseg; // connecting
1249
1250 Int_t c = GetBasicColor();
1251 Int_t i, j;
1252 Int_t indx;
1253 indx = 0;
1254 // outside sphere
1255 // loop all segments on latitudes (except 0 and 180 degrees)
1256 // [0, nlat*fNseg)
1257 Int_t indpar = 0;
1258 for (i = 0; i < nlat; i++) {
1259 for (j = 0; j < fNseg; j++) {
1260 buff.fSegs[indx++] = c;
1261 buff.fSegs[indx++] = i * nlong + j;
1262 buff.fSegs[indx++] = i * nlong + (j + 1) % nlong;
1263 }
1264 }
1265 // loop all segments on longitudes
1266 // nlat*fNseg + [0, (nlat-1)*nlong)
1268 for (i = 0; i < nlat - 1; i++) {
1269 for (j = 0; j < nlong; j++) {
1270 buff.fSegs[indx++] = c;
1271 buff.fSegs[indx++] = i * nlong + j;
1272 buff.fSegs[indx++] = (i + 1) * nlong + j;
1273 }
1274 }
1275 Int_t indup = indlong + (nlat - 1) * nlong;
1276 // extra longitudes on top
1277 // nlat*fNseg+(nlat-1)*nlong + [0, nlong)
1278 if (nup) {
1279 Int_t indpup = nlat * nlong;
1280 for (j = 0; j < nlong; j++) {
1281 buff.fSegs[indx++] = c;
1282 buff.fSegs[indx++] = j;
1283 buff.fSegs[indx++] = indpup;
1284 }
1285 }
1286 Int_t inddown = indup + nup * nlong;
1287 // extra longitudes on bottom
1288 // nlat*fNseg+(nlat+nup-1)*nlong + [0, nlong)
1289 if (ndown) {
1290 Int_t indpdown = nlat * nlong + nup;
1291 for (j = 0; j < nlong; j++) {
1292 buff.fSegs[indx++] = c;
1293 buff.fSegs[indx++] = (nlat - 1) * nlong + j;
1294 buff.fSegs[indx++] = indpdown;
1295 }
1296 }
1302 // inner sphere
1303 Int_t indptin = nlat * nlong + nup + ndown;
1305 // nlat*fNseg+(nlat+nup+ndown-1)*nlong
1306 if (TestShapeBit(kGeoRSeg)) {
1308 indupin = indlongin + (nlat - 1) * nlong;
1309 inddownin = indupin + nup * nlong;
1310 // loop all segments on latitudes (except 0 and 180 degrees)
1311 // indsegin + [0, nlat*fNseg)
1312 for (i = 0; i < nlat; i++) {
1313 for (j = 0; j < fNseg; j++) {
1314 buff.fSegs[indx++] = c + 1;
1315 buff.fSegs[indx++] = indptin + i * nlong + j;
1316 buff.fSegs[indx++] = indptin + i * nlong + (j + 1) % nlong;
1317 }
1318 }
1319 // loop all segments on longitudes
1320 // indsegin + nlat*fNseg + [0, (nlat-1)*nlong)
1321 for (i = 0; i < nlat - 1; i++) {
1322 for (j = 0; j < nlong; j++) {
1323 buff.fSegs[indx++] = c + 1;
1324 buff.fSegs[indx++] = indptin + i * nlong + j;
1325 buff.fSegs[indx++] = indptin + (i + 1) * nlong + j;
1326 }
1327 }
1328 // extra longitudes on top
1329 // indsegin + nlat*fNseg+(nlat-1)*nlong + [0, nlong)
1330 if (nup) {
1332 for (j = 0; j < nlong; j++) {
1333 buff.fSegs[indx++] = c + 1;
1334 buff.fSegs[indx++] = indptin + j;
1335 buff.fSegs[indx++] = indupltop;
1336 }
1337 }
1338 // extra longitudes on bottom
1339 // indsegin + nlat*fNseg+(nlat+nup-1)*nlong + [0, nlong)
1340 if (ndown) {
1342 for (j = 0; j < nlong; j++) {
1343 buff.fSegs[indx++] = c + 1;
1344 buff.fSegs[indx++] = indptin + (nlat - 1) * nlong + j;
1345 buff.fSegs[indx++] = indpdown;
1346 }
1347 }
1349 }
1351 // Segments on phi planes
1352 if (TestShapeBit(kGeoPhiSeg)) {
1353 indtheta += 2 * nlat + nup + ndown;
1354 for (j = 0; j < nlat; j++) {
1355 buff.fSegs[indx++] = c + 2;
1356 buff.fSegs[indx++] = j * nlong;
1358 buff.fSegs[indx++] = indptin + j * nlong;
1359 else
1360 buff.fSegs[indx++] = iptcenter;
1361 }
1362 for (j = 0; j < nlat; j++) {
1363 buff.fSegs[indx++] = c + 2;
1364 buff.fSegs[indx++] = (j + 1) * nlong - 1;
1366 buff.fSegs[indx++] = indptin + (j + 1) * nlong - 1;
1367 else
1368 buff.fSegs[indx++] = iptcenter;
1369 }
1370 if (nup) {
1371 buff.fSegs[indx++] = c + 2;
1372 buff.fSegs[indx++] = nlat * nlong;
1374 buff.fSegs[indx++] = indptin + nlat * nlong;
1375 else
1376 buff.fSegs[indx++] = iptcenter;
1377 }
1378 if (ndown) {
1379 buff.fSegs[indx++] = c + 2;
1380 buff.fSegs[indx++] = nlat * nlong + nup;
1382 buff.fSegs[indx++] = indptin + nlat * nlong + nup;
1383 else
1384 buff.fSegs[indx++] = iptcenter;
1385 }
1386 }
1387 // Segments on cones
1388 if (!nup) {
1389 for (j = 0; j < nlong; j++) {
1390 buff.fSegs[indx++] = c + 2;
1391 buff.fSegs[indx++] = j;
1393 buff.fSegs[indx++] = indptin + j;
1394 else
1395 buff.fSegs[indx++] = iptcenter;
1396 }
1397 }
1398 if (!ndown) {
1399 for (j = 0; j < nlong; j++) {
1400 buff.fSegs[indx++] = c + 2;
1401 buff.fSegs[indx++] = (nlat - 1) * nlong + j;
1403 buff.fSegs[indx++] = indptin + (nlat - 1) * nlong + j;
1404 else
1405 buff.fSegs[indx++] = iptcenter;
1406 }
1407 }
1408
1409 indx = 0;
1410 // Fill polygons for outside sphere (except 0/180)
1411 for (i = 0; i < nlat - 1; i++) {
1412 for (j = 0; j < fNseg; j++) {
1413 buff.fPols[indx++] = c;
1414 buff.fPols[indx++] = 4;
1415 buff.fPols[indx++] = indpar + i * fNseg + j;
1416 buff.fPols[indx++] = indlong + i * nlong + (j + 1) % nlong;
1417 buff.fPols[indx++] = indpar + (i + 1) * fNseg + j;
1418 buff.fPols[indx++] = indlong + i * nlong + j;
1419 }
1420 }
1421 // upper
1422 if (nup) {
1423 for (j = 0; j < fNseg; j++) {
1424 buff.fPols[indx++] = c;
1425 buff.fPols[indx++] = 3;
1426 buff.fPols[indx++] = indup + j;
1427 buff.fPols[indx++] = indup + (j + 1) % nlong;
1428 buff.fPols[indx++] = indpar + j;
1429 }
1430 }
1431 // lower
1432 if (ndown) {
1433 for (j = 0; j < fNseg; j++) {
1434 buff.fPols[indx++] = c;
1435 buff.fPols[indx++] = 3;
1436 buff.fPols[indx++] = inddown + j;
1437 buff.fPols[indx++] = indpar + (nlat - 1) * fNseg + j;
1438 buff.fPols[indx++] = inddown + (j + 1) % nlong;
1439 }
1440 }
1441 // Fill polygons for inside sphere (except 0/180)
1442
1443 if (TestShapeBit(kGeoRSeg)) {
1444 for (i = 0; i < nlat - 1; i++) {
1445 for (j = 0; j < fNseg; j++) {
1446 buff.fPols[indx++] = c + 1;
1447 buff.fPols[indx++] = 4;
1448 buff.fPols[indx++] = indparin + i * fNseg + j;
1449 buff.fPols[indx++] = indlongin + i * nlong + j;
1450 buff.fPols[indx++] = indparin + (i + 1) * fNseg + j;
1451 buff.fPols[indx++] = indlongin + i * nlong + (j + 1) % nlong;
1452 }
1453 }
1454 // upper
1455 if (nup) {
1456 for (j = 0; j < fNseg; j++) {
1457 buff.fPols[indx++] = c + 1;
1458 buff.fPols[indx++] = 3;
1459 buff.fPols[indx++] = indupin + j;
1460 buff.fPols[indx++] = indparin + j;
1461 buff.fPols[indx++] = indupin + (j + 1) % nlong;
1462 }
1463 }
1464 // lower
1465 if (ndown) {
1466 for (j = 0; j < fNseg; j++) {
1467 buff.fPols[indx++] = c + 1;
1468 buff.fPols[indx++] = 3;
1469 buff.fPols[indx++] = inddownin + j;
1470 buff.fPols[indx++] = inddownin + (j + 1) % nlong;
1471 buff.fPols[indx++] = indparin + (nlat - 1) * fNseg + j;
1472 }
1473 }
1474 }
1475 // Polygons on phi planes
1476 if (TestShapeBit(kGeoPhiSeg)) {
1477 for (i = 0; i < nlat - 1; i++) {
1478 buff.fPols[indx++] = c + 2;
1479 if (TestShapeBit(kGeoRSeg)) {
1480 buff.fPols[indx++] = 4;
1481 buff.fPols[indx++] = indlong + i * nlong;
1482 buff.fPols[indx++] = indphi + i + 1;
1483 buff.fPols[indx++] = indlongin + i * nlong;
1484 buff.fPols[indx++] = indphi + i;
1485 } else {
1486 buff.fPols[indx++] = 3;
1487 buff.fPols[indx++] = indlong + i * nlong;
1488 buff.fPols[indx++] = indphi + i + 1;
1489 buff.fPols[indx++] = indphi + i;
1490 }
1491 }
1492 for (i = 0; i < nlat - 1; i++) {
1493 buff.fPols[indx++] = c + 2;
1494 if (TestShapeBit(kGeoRSeg)) {
1495 buff.fPols[indx++] = 4;
1496 buff.fPols[indx++] = indlong + (i + 1) * nlong - 1;
1497 buff.fPols[indx++] = indphi + nlat + i;
1498 buff.fPols[indx++] = indlongin + (i + 1) * nlong - 1;
1499 buff.fPols[indx++] = indphi + nlat + i + 1;
1500 } else {
1501 buff.fPols[indx++] = 3;
1502 buff.fPols[indx++] = indlong + (i + 1) * nlong - 1;
1503 buff.fPols[indx++] = indphi + nlat + i;
1504 buff.fPols[indx++] = indphi + nlat + i + 1;
1505 }
1506 }
1507 if (nup) {
1508 buff.fPols[indx++] = c + 2;
1509 if (TestShapeBit(kGeoRSeg)) {
1510 buff.fPols[indx++] = 4;
1511 buff.fPols[indx++] = indup;
1512 buff.fPols[indx++] = indphi;
1513 buff.fPols[indx++] = indupin;
1514 buff.fPols[indx++] = indphi + 2 * nlat;
1515 } else {
1516 buff.fPols[indx++] = 3;
1517 buff.fPols[indx++] = indup;
1518 buff.fPols[indx++] = indphi;
1519 buff.fPols[indx++] = indphi + 2 * nlat;
1520 }
1521 buff.fPols[indx++] = c + 2;
1522 if (TestShapeBit(kGeoRSeg)) {
1523 buff.fPols[indx++] = 4;
1524 buff.fPols[indx++] = indup + nlong - 1;
1525 buff.fPols[indx++] = indphi + 2 * nlat;
1526 buff.fPols[indx++] = indupin + nlong - 1;
1527 buff.fPols[indx++] = indphi + nlat;
1528 } else {
1529 buff.fPols[indx++] = 3;
1530 buff.fPols[indx++] = indup + nlong - 1;
1531 buff.fPols[indx++] = indphi + 2 * nlat;
1532 buff.fPols[indx++] = indphi + nlat;
1533 }
1534 }
1535 if (ndown) {
1536 buff.fPols[indx++] = c + 2;
1537 if (TestShapeBit(kGeoRSeg)) {
1538 buff.fPols[indx++] = 4;
1539 buff.fPols[indx++] = inddown;
1540 buff.fPols[indx++] = indphi + 2 * nlat + nup;
1541 buff.fPols[indx++] = inddownin;
1542 buff.fPols[indx++] = indphi + nlat - 1;
1543 } else {
1544 buff.fPols[indx++] = 3;
1545 buff.fPols[indx++] = inddown;
1546 buff.fPols[indx++] = indphi + 2 * nlat + nup;
1547 buff.fPols[indx++] = indphi + nlat - 1;
1548 }
1549 buff.fPols[indx++] = c + 2;
1550 if (TestShapeBit(kGeoRSeg)) {
1551 buff.fPols[indx++] = 4;
1552 buff.fPols[indx++] = inddown + nlong - 1;
1553 buff.fPols[indx++] = indphi + 2 * nlat - 1;
1554 buff.fPols[indx++] = inddownin + nlong - 1;
1555 buff.fPols[indx++] = indphi + 2 * nlat + nup;
1556 } else {
1557 buff.fPols[indx++] = 3;
1558 buff.fPols[indx++] = inddown + nlong - 1;
1559 buff.fPols[indx++] = indphi + 2 * nlat - 1;
1560 buff.fPols[indx++] = indphi + 2 * nlat + nup;
1561 }
1562 }
1563 }
1564 // Polygons on cones
1565 if (!nup) {
1566 for (j = 0; j < fNseg; j++) {
1567 buff.fPols[indx++] = c + 2;
1568 if (TestShapeBit(kGeoRSeg)) {
1569 buff.fPols[indx++] = 4;
1570 buff.fPols[indx++] = indpar + j;
1571 buff.fPols[indx++] = indtheta + j;
1572 buff.fPols[indx++] = indparin + j;
1573 buff.fPols[indx++] = indtheta + (j + 1) % nlong;
1574 } else {
1575 buff.fPols[indx++] = 3;
1576 buff.fPols[indx++] = indpar + j;
1577 buff.fPols[indx++] = indtheta + j;
1578 buff.fPols[indx++] = indtheta + (j + 1) % nlong;
1579 }
1580 }
1581 }
1582 if (!ndown) {
1583 for (j = 0; j < fNseg; j++) {
1584 buff.fPols[indx++] = c + 2;
1585 if (TestShapeBit(kGeoRSeg)) {
1586 buff.fPols[indx++] = 4;
1587 buff.fPols[indx++] = indpar + (nlat - 1) * fNseg + j;
1588 buff.fPols[indx++] = indtheta + (1 - nup) * nlong + (j + 1) % nlong;
1589 buff.fPols[indx++] = indparin + (nlat - 1) * fNseg + j;
1590 buff.fPols[indx++] = indtheta + (1 - nup) * nlong + j;
1591 } else {
1592 buff.fPols[indx++] = 3;
1593 buff.fPols[indx++] = indpar + (nlat - 1) * fNseg + j;
1594 buff.fPols[indx++] = indtheta + (1 - nup) * nlong + (j + 1) % nlong;
1595 buff.fPols[indx++] = indtheta + (1 - nup) * nlong + j;
1596 }
1597 }
1598 }
1599}
1600
1601////////////////////////////////////////////////////////////////////////////////
1602/// computes the closest distance from given point to this shape, according
1603/// to option. The matching point on the shape is stored in spoint.
1604
1606{
1607 Double_t r2 = point[0] * point[0] + point[1] * point[1] + point[2] * point[2];
1610 if (r <= 1E-20)
1611 rzero = kTRUE;
1612 // localize theta
1613 Double_t th = 0.;
1614 if (TestShapeBit(kGeoThetaSeg) && (!rzero)) {
1615 th = TMath::ACos(point[2] / r) * TMath::RadToDeg();
1616 }
1617 Double_t saf[4];
1619 ? TGeoShape::Big()
1620 : r - fRmin;
1621 saf[1] = fRmax - r;
1622 saf[2] = saf[3] = TGeoShape::Big();
1624 if (fTheta1 > 0)
1625 saf[2] = r * TMath::Sin((th - fTheta1) * TMath::DegToRad());
1626 if (fTheta2 < 180)
1627 saf[3] = r * TMath::Sin((fTheta2 - th) * TMath::DegToRad());
1628 }
1631 safphi = TGeoShape::SafetyPhi(point, in, fPhi1, fPhi2);
1632 if (in) {
1634 return TMath::Min(safe, safphi);
1635 }
1636 for (Int_t i = 0; i < 4; i++)
1637 saf[i] = -saf[i];
1640 return TMath::Max(safe, safphi);
1641 return safe;
1642}
1643
1644////////////////////////////////////////////////////////////////////////////////
1645/// Save a primitive as a C++ statement(s) on output stream "out".
1646
1647void TGeoSphere::SavePrimitive(std::ostream &out, Option_t * /*option*/ /*= ""*/)
1648{
1650 return;
1651 out << " // Shape: " << GetName() << " type: " << ClassName() << std::endl;
1652 out << " rmin = " << fRmin << ";" << std::endl;
1653 out << " rmax = " << fRmax << ";" << std::endl;
1654 out << " theta1 = " << fTheta1 << ";" << std::endl;
1655 out << " theta2 = " << fTheta2 << ";" << std::endl;
1656 out << " phi1 = " << fPhi1 << ";" << std::endl;
1657 out << " phi2 = " << fPhi2 << ";" << std::endl;
1658 out << " TGeoShape *" << GetPointerName() << " = new TGeoSphere(\"" << GetName()
1659 << "\",rmin,rmax,theta1, theta2,phi1,phi2);" << std::endl;
1661}
1662
1663////////////////////////////////////////////////////////////////////////////////
1664/// Set spherical segment dimensions.
1665
1667 Double_t phi2)
1668{
1669 if (rmin >= rmax) {
1670 Error("SetDimensions", "invalid parameters rmin/rmax");
1671 return;
1672 }
1673 fRmin = rmin;
1674 fRmax = rmax;
1675 if (rmin > 0)
1677 if (theta1 >= theta2 || theta1 < 0 || theta1 > 180 || theta2 > 180) {
1678 Error("SetDimensions", "invalid parameters theta1/theta2");
1679 return;
1680 }
1681 fTheta1 = theta1;
1682 fTheta2 = theta2;
1683 if ((theta2 - theta1) < 180.)
1685 fPhi1 = phi1;
1686 if (phi1 < 0)
1687 fPhi1 += 360.;
1688 fPhi2 = phi2;
1689 while (fPhi2 <= fPhi1)
1690 fPhi2 += 360.;
1693}
1694
1695////////////////////////////////////////////////////////////////////////////////
1696/// Set dimensions of the spherical segment starting from a list of parameters.
1697
1699{
1700 Double_t rmin = param[0];
1701 Double_t rmax = param[1];
1702 Double_t theta1 = 0;
1703 Double_t theta2 = 180.;
1704 Double_t phi1 = 0;
1705 Double_t phi2 = 360.;
1706 if (nparam > 2)
1707 theta1 = param[2];
1708 if (nparam > 3)
1709 theta2 = param[3];
1710 if (nparam > 4)
1711 phi1 = param[4];
1712 if (nparam > 5)
1713 phi2 = param[5];
1715}
1716
1717////////////////////////////////////////////////////////////////////////////////
1718/// Set dimensions of the spherical segment starting from a list of parameters.
1719/// Only takes rmin and rmax
1720
1722{
1723 SetDimensions(param, 2);
1724}
1725
1726////////////////////////////////////////////////////////////////////////////////
1727/// Set the number of divisions of mesh circles keeping aspect ratio.
1728
1730{
1731 fNseg = p;
1733 if (dphi < 0)
1734 dphi += 360;
1736 fNz = Int_t(fNseg * dtheta / dphi) + 1;
1737 if (fNz < 2)
1738 fNz = 2;
1739}
1740
1741////////////////////////////////////////////////////////////////////////////////
1742/// create sphere mesh points
1743
1745{
1746 if (!points) {
1747 Error("SetPoints", "Input array is NULL");
1748 return;
1749 }
1750 Bool_t full = kTRUE;
1752 full = kFALSE;
1753 Int_t ncenter = 1;
1754 if (full || TestShapeBit(kGeoRSeg))
1755 ncenter = 0;
1756 Int_t nup = (fTheta1 > 0) ? 0 : 1;
1757 Int_t ndown = (fTheta2 < 180) ? 0 : 1;
1758 // number of different latitudes, excluding 0 and 180 degrees
1759 Int_t nlat = fNz + 1 - (nup + ndown);
1760 // number of different longitudes
1761 Int_t nlong = fNseg;
1763 nlong++;
1764 // total number of points on mesh is:
1765 // nlat*nlong + nup + ndown + ncenter; // in case rmin=0
1766 // 2*(nlat*nlong + nup + ndown); // in case rmin>0
1767 Int_t i, j;
1770 Double_t dphi = (phi2 - phi1) / fNseg;
1774 Double_t z, zi, theta, phi, cphi, sphi;
1775 Int_t indx = 0;
1776 // FILL ALL POINTS ON OUTER SPHERE
1777 // (nlat * nlong) points
1778 // loop all latitudes except 0/180 degrees (nlat times)
1779 // ilat = [0,nlat] jlong = [0,nlong]
1780 // Index(ilat, jlong) = 3*(ilat*nlat + jlong)
1781 for (i = 0; i < nlat; i++) {
1782 theta = theta1 + (nup + i) * dtheta;
1783 z = fRmax * TMath::Cos(theta);
1784 zi = fRmax * TMath::Sin(theta);
1785 // loop all different longitudes (nlong times)
1786 for (j = 0; j < nlong; j++) {
1787 phi = phi1 + j * dphi;
1788 cphi = TMath::Cos(phi);
1789 sphi = TMath::Sin(phi);
1790 points[indx++] = zi * cphi;
1791 points[indx++] = zi * sphi;
1792 points[indx++] = z;
1793 }
1794 }
1795 // upper/lower points (if they exist) for outer sphere
1796 if (nup) {
1797 // ind_up = 3*nlat*nlong
1798 points[indx++] = 0.;
1799 points[indx++] = 0.;
1800 points[indx++] = fRmax;
1801 }
1802 if (ndown) {
1803 // ind_down = 3*(nlat*nlong+nup)
1804 points[indx++] = 0.;
1805 points[indx++] = 0.;
1806 points[indx++] = -fRmax;
1807 }
1808 // do the same for inner sphere if it exist
1809 // Start_index = 3*(nlat*nlong + nup + ndown)
1810 if (TestShapeBit(kGeoRSeg)) {
1811 // Index(ilat, jlong) = start_index + 3*(ilat*nlat + jlong)
1812 for (i = 0; i < nlat; i++) {
1813 theta = theta1 + (nup + i) * dtheta;
1814 z = fRmin * TMath::Cos(theta);
1815 zi = fRmin * TMath::Sin(theta);
1816 // loop all different longitudes (nlong times)
1817 for (j = 0; j < nlong; j++) {
1818 phi = phi1 + j * dphi;
1819 cphi = TMath::Cos(phi);
1820 sphi = TMath::Sin(phi);
1821 points[indx++] = zi * cphi;
1822 points[indx++] = zi * sphi;
1823 points[indx++] = z;
1824 }
1825 }
1826 // upper/lower points (if they exist) for inner sphere
1827 if (nup) {
1828 // ind_up = start_index + 3*nlat*nlong
1829 points[indx++] = 0.;
1830 points[indx++] = 0.;
1831 points[indx++] = fRmin;
1832 }
1833 if (ndown) {
1834 // ind_down = start_index + 3*(nlat*nlong+nup)
1835 points[indx++] = 0.;
1836 points[indx++] = 0.;
1837 points[indx++] = -fRmin;
1838 }
1839 }
1840 // Add center of sphere if needed
1841 if (ncenter) {
1842 // ind_center = 6*(nlat*nlong + nup + ndown)
1843 points[indx++] = 0.;
1844 points[indx++] = 0.;
1845 points[indx++] = 0.;
1846 }
1847}
1848
1849////////////////////////////////////////////////////////////////////////////////
1850/// create sphere mesh points
1851
1853{
1854 if (!points) {
1855 Error("SetPoints", "Input array is NULL");
1856 return;
1857 }
1858 Bool_t full = kTRUE;
1860 full = kFALSE;
1861 Int_t ncenter = 1;
1862 if (full || TestShapeBit(kGeoRSeg))
1863 ncenter = 0;
1864 Int_t nup = (fTheta1 > 0) ? 0 : 1;
1865 Int_t ndown = (fTheta2 < 180) ? 0 : 1;
1866 // number of different latitudes, excluding 0 and 180 degrees
1867 Int_t nlat = fNz + 1 - (nup + ndown);
1868 // number of different longitudes
1869 Int_t nlong = fNseg;
1871 nlong++;
1872 // total number of points on mesh is:
1873 // nlat*nlong + nup + ndown + ncenter; // in case rmin=0
1874 // 2*(nlat*nlong + nup + ndown); // in case rmin>0
1875 Int_t i, j;
1878 Double_t dphi = (phi2 - phi1) / fNseg;
1882 Double_t z, zi, theta, phi, cphi, sphi;
1883 Int_t indx = 0;
1884 // FILL ALL POINTS ON OUTER SPHERE
1885 // (nlat * nlong) points
1886 // loop all latitudes except 0/180 degrees (nlat times)
1887 // ilat = [0,nlat] jlong = [0,nlong]
1888 // Index(ilat, jlong) = 3*(ilat*nlat + jlong)
1889 for (i = 0; i < nlat; i++) {
1890 theta = theta1 + (nup + i) * dtheta;
1891 z = fRmax * TMath::Cos(theta);
1892 zi = fRmax * TMath::Sin(theta);
1893 // loop all different longitudes (nlong times)
1894 for (j = 0; j < nlong; j++) {
1895 phi = phi1 + j * dphi;
1896 cphi = TMath::Cos(phi);
1897 sphi = TMath::Sin(phi);
1898 points[indx++] = zi * cphi;
1899 points[indx++] = zi * sphi;
1900 points[indx++] = z;
1901 }
1902 }
1903 // upper/lower points (if they exist) for outer sphere
1904 if (nup) {
1905 // ind_up = 3*nlat*nlong
1906 points[indx++] = 0.;
1907 points[indx++] = 0.;
1908 points[indx++] = fRmax;
1909 }
1910 if (ndown) {
1911 // ind_down = 3*(nlat*nlong+nup)
1912 points[indx++] = 0.;
1913 points[indx++] = 0.;
1914 points[indx++] = -fRmax;
1915 }
1916 // do the same for inner sphere if it exist
1917 // Start_index = 3*(nlat*nlong + nup + ndown)
1918 if (TestShapeBit(kGeoRSeg)) {
1919 // Index(ilat, jlong) = start_index + 3*(ilat*nlat + jlong)
1920 for (i = 0; i < nlat; i++) {
1921 theta = theta1 + (nup + i) * dtheta;
1922 z = fRmin * TMath::Cos(theta);
1923 zi = fRmin * TMath::Sin(theta);
1924 // loop all different longitudes (nlong times)
1925 for (j = 0; j < nlong; j++) {
1926 phi = phi1 + j * dphi;
1927 cphi = TMath::Cos(phi);
1928 sphi = TMath::Sin(phi);
1929 points[indx++] = zi * cphi;
1930 points[indx++] = zi * sphi;
1931 points[indx++] = z;
1932 }
1933 }
1934 // upper/lower points (if they exist) for inner sphere
1935 if (nup) {
1936 // ind_up = start_index + 3*nlat*nlong
1937 points[indx++] = 0.;
1938 points[indx++] = 0.;
1939 points[indx++] = fRmin;
1940 }
1941 if (ndown) {
1942 // ind_down = start_index + 3*(nlat*nlong+nup)
1943 points[indx++] = 0.;
1944 points[indx++] = 0.;
1945 points[indx++] = -fRmin;
1946 }
1947 }
1948 // Add center of sphere if needed
1949 if (ncenter) {
1950 // ind_center = 6*(nlat*nlong + nup + ndown)
1951 points[indx++] = 0.;
1952 points[indx++] = 0.;
1953 points[indx++] = 0.;
1954 }
1955}
1956
1957////////////////////////////////////////////////////////////////////////////////
1958/// Returns numbers of vertices, segments and polygons composing the shape mesh.
1959
1961{
1962 TGeoSphere *localThis = const_cast<TGeoSphere *>(this);
1963 localThis->SetNumberOfDivisions(gGeoManager->GetNsegments());
1964 Bool_t full = kTRUE;
1966 full = kFALSE;
1967 Int_t ncenter = 1;
1968 if (full || TestShapeBit(kGeoRSeg))
1969 ncenter = 0;
1970 Int_t nup = (fTheta1 > 0) ? 0 : 1;
1971 Int_t ndown = (fTheta2 < 180) ? 0 : 1;
1972 // number of different latitudes, excluding 0 and 180 degrees
1973 Int_t nlat = fNz + 1 - (nup + ndown);
1974 // number of different longitudes
1975 Int_t nlong = fNseg;
1977 nlong++;
1978
1979 nvert = nlat * nlong + nup + ndown + ncenter;
1981 nvert *= 2;
1982
1983 nsegs = nlat * fNseg + (nlat - 1 + nup + ndown) * nlong; // outer sphere
1985 nsegs *= 2; // inner sphere
1987 nsegs += 2 * nlat + nup + ndown; // 2 phi planes
1988 nsegs += nlong * (2 - nup - ndown); // connecting cones
1989
1990 npols = fNz * fNseg; // outer
1992 npols *= 2; // inner
1994 npols += 2 * fNz; // 2 phi planes
1995 npols += (2 - nup - ndown) * fNseg; // connecting
1996}
1997
1998////////////////////////////////////////////////////////////////////////////////
1999/// Return number of vertices of the mesh representation
2000
2002{
2003 Bool_t full = kTRUE;
2005 full = kFALSE;
2006 Int_t ncenter = 1;
2007 if (full || TestShapeBit(kGeoRSeg))
2008 ncenter = 0;
2009 Int_t nup = (fTheta1 > 0) ? 0 : 1;
2010 Int_t ndown = (fTheta2 < 180) ? 0 : 1;
2011 // number of different latitudes, excluding 0 and 180 degrees
2012 Int_t nlat = fNz + 1 - (nup + ndown);
2013 // number of different longitudes
2014 Int_t nlong = fNseg;
2016 nlong++;
2017 // total number of points on mesh is:
2018 // nlat*nlong + nup + ndown + ncenter; // in case rmin=0
2019 // 2*(nlat*nlong + nup + ndown); // in case rmin>0
2020 Int_t numPoints = 0;
2022 numPoints = 2 * (nlat * nlong + nup + ndown);
2023 else
2024 numPoints = nlat * nlong + nup + ndown + ncenter;
2025 return numPoints;
2026}
2027
2028////////////////////////////////////////////////////////////////////////////////
2029////// obsolete - to be removed
2030
2032
2033////////////////////////////////////////////////////////////////////////////////
2034/// Fills a static 3D buffer and returns a reference.
2035
2037{
2038 static TBuffer3DSphere buffer;
2039
2041
2043 buffer.fRadiusInner = fRmin;
2044 buffer.fRadiusOuter = fRmax;
2045 buffer.fThetaMin = fTheta1;
2046 buffer.fThetaMax = fTheta2;
2047 buffer.fPhiMin = fPhi1;
2048 buffer.fPhiMax = fPhi2;
2050 }
2052 // We want FillBuffer to be const
2053 TGeoSphere *localThis = const_cast<TGeoSphere *>(this);
2054 localThis->SetNumberOfDivisions(gGeoManager->GetNsegments());
2055
2056 Bool_t full = kTRUE;
2058 full = kFALSE;
2059 Int_t ncenter = 1;
2060 if (full || TestShapeBit(kGeoRSeg))
2061 ncenter = 0;
2062 Int_t nup = (fTheta1 > 0) ? 0 : 1;
2063 Int_t ndown = (fTheta2 < 180) ? 0 : 1;
2064 // number of different latitudes, excluding 0 and 180 degrees
2065 Int_t nlat = fNz + 1 - (nup + ndown);
2066 // number of different longitudes
2067 Int_t nlong = fNseg;
2069 nlong++;
2070
2073 nbPnts *= 2;
2074
2075 Int_t nbSegs = nlat * fNseg + (nlat - 1 + nup + ndown) * nlong; // outer sphere
2077 nbSegs *= 2; // inner sphere
2079 nbSegs += 2 * nlat + nup + ndown; // 2 phi planes
2080 nbSegs += nlong * (2 - nup - ndown); // connecting cones
2081
2082 Int_t nbPols = fNz * fNseg; // outer
2084 nbPols *= 2; // inner
2086 nbPols += 2 * fNz; // 2 phi planes
2087 nbPols += (2 - nup - ndown) * fNseg; // connecting
2088
2089 if (buffer.SetRawSizes(nbPnts, 3 * nbPnts, nbSegs, 3 * nbSegs, nbPols, 6 * nbPols)) {
2091 }
2092 }
2094 SetPoints(buffer.fPnts);
2095 if (!buffer.fLocalFrame) {
2096 TransformPoints(buffer.fPnts, buffer.NbPnts());
2097 }
2098 SetSegsAndPols(buffer);
2100 }
2101
2102 return buffer;
2103}
2104
2105////////////////////////////////////////////////////////////////////////////////
2106/// Check the inside status for each of the points in the array.
2107/// Input: Array of point coordinates + vector size
2108/// Output: Array of Booleans for the inside of each point
2109
2111{
2112 for (Int_t i = 0; i < vecsize; i++)
2113 inside[i] = Contains(&points[3 * i]);
2114}
2115
2116////////////////////////////////////////////////////////////////////////////////
2117/// Compute the normal for an array o points so that norm.dot.dir is positive
2118/// Input: Arrays of point coordinates and directions + vector size
2119/// Output: Array of normal directions
2120
2122{
2123 for (Int_t i = 0; i < vecsize; i++)
2124 ComputeNormal(&points[3 * i], &dirs[3 * i], &norms[3 * i]);
2125}
2126
2127////////////////////////////////////////////////////////////////////////////////
2128/// Compute distance from array of input points having directions specified by dirs. Store output in dists
2129
2131 Double_t *step) const
2132{
2133 for (Int_t i = 0; i < vecsize; i++)
2134 dists[i] = DistFromInside(&points[3 * i], &dirs[3 * i], 3, step[i]);
2135}
2136
2137////////////////////////////////////////////////////////////////////////////////
2138/// Compute distance from array of input points having directions specified by dirs. Store output in dists
2139
2141 Double_t *step) const
2142{
2143 for (Int_t i = 0; i < vecsize; i++)
2144 dists[i] = DistFromOutside(&points[3 * i], &dirs[3 * i], 3, step[i]);
2145}
2146
2147////////////////////////////////////////////////////////////////////////////////
2148/// Compute safe distance from each of the points in the input array.
2149/// Input: Array of point coordinates, array of statuses for these points, size of the arrays
2150/// Output: Safety values
2151
2153{
2154 for (Int_t i = 0; i < vecsize; i++)
2155 safe[i] = Safety(&points[3 * i], inside[i]);
2156}
#define d(i)
Definition RSha256.hxx:102
#define b(i)
Definition RSha256.hxx:100
#define c(i)
Definition RSha256.hxx:101
#define a(i)
Definition RSha256.hxx:99
#define s1(x)
Definition RSha256.hxx:91
std::size_t capacity
int Int_t
Signed integer 4 bytes (int)
Definition RtypesCore.h:60
float Float_t
Float 4 bytes (float)
Definition RtypesCore.h:72
constexpr Bool_t kFALSE
Definition RtypesCore.h:109
double Double_t
Double 8 bytes.
Definition RtypesCore.h:74
constexpr Bool_t kTRUE
Definition RtypesCore.h:108
const char Option_t
Option string (const char)
Definition RtypesCore.h:81
ROOT::Detail::TRangeCast< T, true > TRangeDynCast
TRangeDynCast is an adapter class that allows the typed iteration through a TCollection.
winID h TVirtualViewer3D TVirtualGLPainter p
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t r
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize id
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:142
R__EXTERN TGeoManager * gGeoManager
float xmin
float ymin
float xmax
float ymax
Sphere description class - see TBuffer3DTypes for producer classes Supports hollow and cut spheres.
Definition TBuffer3D.h:139
Double_t fRadiusInner
Definition TBuffer3D.h:153
Double_t fThetaMin
Definition TBuffer3D.h:155
Double_t fPhiMin
Definition TBuffer3D.h:157
Double_t fThetaMax
Definition TBuffer3D.h:156
Double_t fRadiusOuter
Definition TBuffer3D.h:154
Double_t fPhiMax
Definition TBuffer3D.h:158
Generic 3D primitive description class.
Definition TBuffer3D.h:18
UInt_t NbPnts() const
Definition TBuffer3D.h:89
Bool_t SectionsValid(UInt_t mask) const
Definition TBuffer3D.h:76
@ kShapeSpecific
Definition TBuffer3D.h:61
void SetSectionsValid(UInt_t mask)
Definition TBuffer3D.h:74
Bool_t fLocalFrame
Definition TBuffer3D.h:99
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:122
Box class.
Definition TGeoBBox.h:18
void FillBuffer3D(TBuffer3D &buffer, Int_t reqSections, Bool_t localFrame) const override
Fills the supplied buffer, with sections in desired frame See TBuffer3D.h for explanation of sections...
Double_t fDX
Definition TGeoBBox.h:21
void SetBoxDimensions(Double_t dx, Double_t dy, Double_t dz, Double_t *origin=nullptr)
Set parameters of the box.
Double_t DistFromOutside(const Double_t *point, const Double_t *dir, Int_t iact=1, Double_t step=TGeoShape::Big(), Double_t *safe=nullptr) const override
Compute distance from outside point to surface of the box.
Definition TGeoBBox.cxx:431
Double_t fOrigin[3]
Definition TGeoBBox.h:24
void InspectShape() const override
Prints shape parameters.
Definition TGeoBBox.cxx:845
Double_t fDY
Definition TGeoBBox.h:22
Double_t fDZ
Definition TGeoBBox.h:23
static void DistToCone(const Double_t *point, const Double_t *dir, Double_t dz, Double_t r1, Double_t r2, Double_t &b, Double_t &delta)
Static method to compute distance to a conical surface with :
Definition TGeoCone.cxx:564
TGeoVolumeMulti * MakeVolumeMulti(const char *name, TGeoMedium *medium)
Make a TGeoVolumeMulti handling a list of volumes.
Int_t GetNsegments() const
Get number of segments approximating circles.
Node containing an offset.
Definition TGeoNode.h:185
base finder class for patterns. A pattern is specifying a division type
a spherical phi divison pattern
a spherical R divison pattern
a spherical theta divison pattern
Base abstract class for all shapes.
Definition TGeoShape.h:25
static Double_t Big()
Definition TGeoShape.h:95
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 DistToPhiMin(const Double_t *point, const Double_t *dir, Double_t s1, Double_t c1, Double_t s2, Double_t c2, Double_t sm, Double_t cm, Bool_t in=kTRUE)
compute distance from point (inside phi) to both phi planes. Return minimum.
static Double_t SafetyPhi(const Double_t *point, Bool_t in, Double_t phi1, Double_t phi2)
Static method to compute safety w.r.t a phi corner defined by cosines/sines of the angles phi1,...
static Bool_t IsSameWithinTolerance(Double_t a, Double_t b)
Check if two numbers differ with less than a tolerance.
const char * GetPointerName() const
Provide a pointer name containing uid.
Int_t ShapeDistancetoPrimitive(Int_t numpoints, Int_t px, Int_t py) const
Returns distance to shape primitive mesh.
static void NormalPhi(const Double_t *point, const Double_t *dir, Double_t *norm, Double_t c1, Double_t s1, Double_t c2, Double_t s2)
Static method to compute normal to phi planes.
const char * GetName() const override
Get the shape name.
@ kGeoSavePrimitive
Definition TGeoShape.h:65
@ kGeoThetaSeg
Definition TGeoShape.h:36
static Double_t Tolerance()
Definition TGeoShape.h:98
static Bool_t IsCloseToPhi(Double_t epsil, const Double_t *point, Double_t c1, Double_t s1, Double_t c2, Double_t s2)
True if point is closer than epsil to one of the phi planes defined by c1,s1 or c2,...
Bool_t TestShapeBit(UInt_t f) const
Definition TGeoShape.h:177
TGeoSphere are not just balls having internal and external radii, but sectors of a sphere having defi...
Definition TGeoSphere.h:17
Double_t DistToSphere(const Double_t *point, const Double_t *dir, Double_t rsph, Bool_t check=kTRUE, Bool_t firstcross=kTRUE) const
compute distance to sphere of radius rsph. Direction has to be a unit vector
Double_t fPhi1
Definition TGeoSphere.h:25
Double_t fTheta2
Definition TGeoSphere.h:24
TGeoVolume * Divide(TGeoVolume *voldiv, const char *divname, Int_t iaxis, Int_t ndiv, Double_t start, Double_t step) override
Divide this box shape belonging to volume "voldiv" into ndiv equal volumes called divname,...
void DistFromInside_v(const Double_t *points, const Double_t *dirs, Double_t *dists, Int_t vecsize, Double_t *step) const override
Compute distance from array of input points having directions specified by dirs. Store output in dist...
void SetPoints(Double_t *points) const override
create sphere mesh points
Bool_t Contains(const Double_t *point) const override
test if point is inside this sphere check Rmin<=R<=Rmax
void Contains_v(const Double_t *points, Bool_t *inside, Int_t vecsize) const override
Check the inside status for each of the points in the array.
void SetSphDimensions(Double_t rmin, Double_t rmax, Double_t theta1, Double_t theta2, Double_t phi1, Double_t phi2)
Set spherical segment dimensions.
void GetMeshNumbers(Int_t &nvert, Int_t &nsegs, Int_t &npols) const override
Returns numbers of vertices, segments and polygons composing the shape mesh.
Double_t fRmax
Definition TGeoSphere.h:22
void Safety_v(const Double_t *points, const Bool_t *inside, Double_t *safe, Int_t vecsize) const override
Compute safe distance from each of the points in the input array.
Int_t fNseg
Definition TGeoSphere.h:20
void DistFromOutside_v(const Double_t *points, const Double_t *dirs, Double_t *dists, Int_t vecsize, Double_t *step) const override
Compute distance from array of input points having directions specified by dirs. Store output in dist...
Double_t fRmin
Definition TGeoSphere.h:21
TGeoSphere()
Default constructor.
void SetDimensions(Double_t *param) override
Set dimensions of the spherical segment starting from a list of parameters.
TBuffer3D * MakeBuffer3D() const override
Creates a TBuffer3D describing this shape.
void InspectShape() const override
print shape parameters
void ComputeNormal(const Double_t *point, const Double_t *dir, Double_t *norm) const override
Compute normal to closest surface from POINT.
void SetSegsAndPols(TBuffer3D &buff) const override
Fill TBuffer3D structure for segments and polygons.
Bool_t IsPointInside(const Double_t *point, Bool_t checkR=kTRUE, Bool_t checkTh=kTRUE, Bool_t checkPh=kTRUE) const
Check if a point is inside radius/theta/phi ranges for the spherical sector.
Double_t Safety(const Double_t *point, Bool_t in=kTRUE) const override
computes the closest distance from given point to this shape, according to option.
void ComputeBBox() override
compute bounding box of the sphere
Double_t fTheta1
Definition TGeoSphere.h:23
Double_t GetAxisRange(Int_t iaxis, Double_t &xlo, Double_t &xhi) const override
Get range of shape for a given axis.
Double_t DistFromOutside(const Double_t *point, const Double_t *dir, Int_t iact=1, Double_t step=TGeoShape::Big(), Double_t *safe=nullptr) const override
compute distance from outside point to surface of the sphere Check if the bounding box is crossed wit...
const TBuffer3D & GetBuffer3D(Int_t reqSections, Bool_t localFrame) const override
Fills a static 3D buffer and returns a reference.
const char * GetAxisName(Int_t iaxis) const override
Returns name of axis IAXIS.
Double_t Capacity() const override
Computes capacity of the shape in [length^3].
virtual void SetNumberOfDivisions(Int_t p)
Set the number of divisions of mesh circles keeping aspect ratio.
void ComputeNormal_v(const Double_t *points, const Double_t *dirs, Double_t *norms, Int_t vecsize) override
Compute the normal for an array o points so that norm.dot.dir is positive Input: Arrays of point coor...
void GetBoundingCylinder(Double_t *param) const override
Fill vector param[4] with the bounding cylinder parameters.
Int_t GetNmeshVertices() const override
Return number of vertices of the mesh representation.
~TGeoSphere() override
destructor
Double_t DistFromInside(const Double_t *point, const Double_t *dir, Int_t iact=1, Double_t step=TGeoShape::Big(), Double_t *safe=nullptr) const override
compute distance from inside point to surface of the sphere
Int_t IsOnBoundary(const Double_t *point) const
Check if a point in local sphere coordinates is close to a boundary within shape tolerance.
void SavePrimitive(std::ostream &out, Option_t *option="") override
Save a primitive as a C++ statement(s) on output stream "out".
Int_t DistancetoPrimitive(Int_t px, Int_t py) override
compute closest distance from point px,py to each corner
Double_t fPhi2
Definition TGeoSphere.h:26
Int_t fNz
Definition TGeoSphere.h:19
void Sizeof3D() const override
Volume families.
Definition TGeoVolume.h:269
TGeoVolume, TGeoVolumeMulti, TGeoVolumeAssembly are the volume classes.
Definition TGeoVolume.h:45
R__ALWAYS_INLINE Bool_t TestBit(UInt_t f) const
Definition TObject.h:204
virtual const char * ClassName() const
Returns name of class to which the object belongs.
Definition TObject.cxx:226
void SetBit(UInt_t f, Bool_t set)
Set or unset the user status bits as specified in f.
Definition TObject.cxx:886
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition TObject.cxx:1096
Basic string class.
Definition TString.h:137
const char * Data() const
Definition TString.h:385
TPaveText * pt
return c1
Definition legend1.C:41
const Int_t n
Definition legend1.C:16
return c2
Definition legend2.C:14
Double_t ACos(Double_t)
Returns the principal value of the arc cosine of x, expressed in radians.
Definition TMath.h:645
Long64_t LocMin(Long64_t n, const T *a)
Returns index of array with the minimum element.
Definition TMath.h:995
Short_t Max(Short_t a, Short_t b)
Returns the largest of a and b.
Definition TMathBase.h:249
T1 Sign(T1 a, T2 b)
Returns a value with the magnitude of a and the sign of b.
Definition TMathBase.h:174
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:659
Long64_t LocMax(Long64_t n, const T *a)
Returns index of array with the maximum element.
Definition TMath.h:1105
constexpr Double_t DegToRad()
Conversion from degree to radian: .
Definition TMath.h:82
Double_t Sqrt(Double_t x)
Returns the square root of x.
Definition TMath.h:675
Short_t Min(Short_t a, Short_t b)
Returns the smallest of a and b.
Definition TMathBase.h:197
Double_t Cos(Double_t)
Returns the cosine of an angle of x radians.
Definition TMath.h:607
constexpr Double_t Pi()
Definition TMath.h:40
Double_t Sin(Double_t)
Returns the sine of an angle of x radians.
Definition TMath.h:601
constexpr Double_t RadToDeg()
Conversion from radian to degree: .
Definition TMath.h:75
Short_t Abs(Short_t d)
Returns the absolute value of parameter Short_t d.
Definition TMathBase.h:122
auto * th2
Definition textalign.C:18
auto * th1
Definition textalign.C:14