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
66TGeoSphere::TGeoSphere()
67{
68 SetShapeBit(TGeoShape::kGeoSph);
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
82TGeoSphere::TGeoSphere(Double_t rmin, Double_t rmax, Double_t theta1, Double_t theta2, Double_t phi1, Double_t phi2)
83 : TGeoBBox(0, 0, 0)
84{
85 SetShapeBit(TGeoShape::kGeoSph);
86 SetSphDimensions(rmin, rmax, theta1, theta2, phi1, phi2);
87 ComputeBBox();
88 SetNumberOfDivisions(20);
89}
90
91////////////////////////////////////////////////////////////////////////////////
92/// Default constructor specifying minimum and maximum radius
93
94TGeoSphere::TGeoSphere(const char *name, Double_t rmin, Double_t rmax, Double_t theta1, Double_t theta2, Double_t phi1,
95 Double_t phi2)
96 : TGeoBBox(name, 0, 0, 0)
97{
98 SetShapeBit(TGeoShape::kGeoSph);
99 SetSphDimensions(rmin, rmax, theta1, theta2, phi1, phi2);
100 ComputeBBox();
101 SetNumberOfDivisions(20);
102}
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
113TGeoSphere::TGeoSphere(Double_t *param, Int_t nparam) : TGeoBBox(0, 0, 0)
114{
115 SetShapeBit(TGeoShape::kGeoSph);
116 SetDimensions(param, nparam);
117 ComputeBBox();
118 SetNumberOfDivisions(20);
119}
120
121////////////////////////////////////////////////////////////////////////////////
122/// destructor
123
124TGeoSphere::~TGeoSphere() {}
125
126////////////////////////////////////////////////////////////////////////////////
127/// Computes capacity of the shape in [length^3]
128
139
140////////////////////////////////////////////////////////////////////////////////
141/// compute bounding box of the sphere
142
143void TGeoSphere::ComputeBBox()
144{
145 if (TGeoShape::IsSameWithinTolerance(TMath::Abs(fTheta2 - fTheta1), 180)) {
146 if (TGeoShape::IsSameWithinTolerance(TMath::Abs(fPhi2 - fPhi1), 360)) {
147 TGeoBBox::SetBoxDimensions(fRmax, fRmax, fRmax);
148 memset(fOrigin, 0, 3 * sizeof(Double_t));
149 return;
150 }
151 }
152 Double_t st1 = TMath::Sin(fTheta1 * TMath::DegToRad());
153 Double_t st2 = TMath::Sin(fTheta2 * TMath::DegToRad());
154 Double_t r1min, r1max, r2min, r2max, rmin, rmax;
155 r1min = TMath::Min(fRmax * st1, fRmax * st2);
156 r1max = TMath::Max(fRmax * st1, fRmax * st2);
157 r2min = TMath::Min(fRmin * st1, fRmin * st2);
158 r2max = TMath::Max(fRmin * st1, fRmin * st2);
159 if (((fTheta1 <= 90) && (fTheta2 >= 90)) || ((fTheta2 <= 90) && (fTheta1 >= 90))) {
160 r1max = fRmax;
161 r2max = fRmin;
162 }
163 rmin = TMath::Min(r1min, r2min);
164 rmax = TMath::Max(r1max, r2max);
165
166 Double_t xc[4];
167 Double_t yc[4];
168 xc[0] = rmax * TMath::Cos(fPhi1 * TMath::DegToRad());
169 yc[0] = rmax * TMath::Sin(fPhi1 * TMath::DegToRad());
170 xc[1] = rmax * TMath::Cos(fPhi2 * TMath::DegToRad());
171 yc[1] = rmax * TMath::Sin(fPhi2 * TMath::DegToRad());
172 xc[2] = rmin * TMath::Cos(fPhi1 * TMath::DegToRad());
173 yc[2] = rmin * TMath::Sin(fPhi1 * TMath::DegToRad());
174 xc[3] = rmin * TMath::Cos(fPhi2 * TMath::DegToRad());
175 yc[3] = rmin * TMath::Sin(fPhi2 * TMath::DegToRad());
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])];
181 Double_t dp = fPhi2 - fPhi1;
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;
212 xc[0] = fRmax * TMath::Cos(fTheta1 * TMath::DegToRad());
213 xc[1] = fRmax * TMath::Cos(fTheta2 * TMath::DegToRad());
214 xc[2] = fRmin * TMath::Cos(fTheta1 * TMath::DegToRad());
215 xc[3] = fRmin * TMath::Cos(fTheta2 * TMath::DegToRad());
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];
235 Double_t r = TMath::Sqrt(r2);
236 Bool_t rzero = kFALSE;
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];
249 saf[0] = (TGeoShape::IsSameWithinTolerance(fRmin, 0) && !TestShapeBit(kGeoThetaSeg) && !TestShapeBit(kGeoPhiSeg))
250 ? TGeoShape::Big()
251 : TMath::Abs(r - fRmin);
252 saf[1] = TMath::Abs(fRmax - r);
253 saf[2] = saf[3] = TGeoShape::Big();
254 if (TestShapeBit(kGeoThetaSeg)) {
255 if (fTheta1 > 0) {
256 saf[2] = r * TMath::Abs(TMath::Sin(th - fTheta1 * TMath::DegToRad()));
257 }
258 if (fTheta2 < 180) {
259 saf[3] = r * TMath::Abs(TMath::Sin(fTheta2 * TMath::DegToRad() - th));
260 }
261 }
262 Int_t i = TMath::LocMin(4, saf);
263 if (TestShapeBit(kGeoPhiSeg)) {
264 Double_t c1 = TMath::Cos(fPhi1 * TMath::DegToRad());
265 Double_t s1 = TMath::Sin(fPhi1 * TMath::DegToRad());
266 Double_t c2 = TMath::Cos(fPhi2 * TMath::DegToRad());
267 Double_t s2 = TMath::Sin(fPhi2 * TMath::DegToRad());
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
300Int_t TGeoSphere::IsOnBoundary(const Double_t *point) const
301{
302 Int_t icode = 0;
303 Double_t tol = TGeoShape::Tolerance();
304 Double_t r2 = point[0] * point[0] + point[1] * point[1] + point[2] * point[2];
305 Double_t drsqout = r2 - fRmax * fRmax;
306 // Test if point is on fRmax boundary
307 if (TMath::Abs(drsqout) < 2. * fRmax * tol)
308 return 2;
309 Double_t drsqin = r2;
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 }
316 if (TestShapeBit(kGeoPhiSeg)) {
317 Double_t phi = TMath::ATan2(point[1], point[0]);
318 if (phi < 0)
319 phi += 2 * TMath::Pi();
320 Double_t phi1 = fPhi1 * TMath::DegToRad();
321 Double_t phi2 = fPhi2 * TMath::DegToRad();
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 }
329 if (TestShapeBit(kGeoThetaSeg)) {
330 Double_t r = TMath::Sqrt(r2);
331 Double_t theta = TMath::ACos(point[2] / r2);
332 Double_t theta1 = fTheta1 * TMath::DegToRad();
333 Double_t theta2 = fTheta2 * TMath::DegToRad();
334 Double_t ddt;
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
352Bool_t TGeoSphere::IsPointInside(const Double_t *point, Bool_t checkR, Bool_t checkTh, Bool_t checkPh) const
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;
363 if (checkPh && TestShapeBit(kGeoPhiSeg)) {
364 Double_t phi = TMath::ATan2(point[1], point[0]) * TMath::RadToDeg();
365 while (phi < fPhi1)
366 phi += 360.;
367 Double_t dphi = fPhi2 - fPhi1;
368 Double_t ddp = phi - fPhi1;
369 if (ddp > dphi)
370 return kFALSE;
371 }
372 if (checkTh && TestShapeBit(kGeoThetaSeg)) {
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
386Bool_t TGeoSphere::Contains(const Double_t *point) const
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
396 if (TestShapeBit(kGeoPhiSeg)) {
397 Double_t phi = TMath::ATan2(point[1], point[0]) * TMath::RadToDeg();
398 if (phi < 0)
399 phi += 360.;
400 Double_t dphi = fPhi2 - fPhi1;
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 }
409 if (TestShapeBit(kGeoThetaSeg)) {
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
422Int_t TGeoSphere::DistancetoPrimitive(Int_t px, Int_t py)
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
434Double_t
435TGeoSphere::DistFromOutside(const Double_t *point, const Double_t *dir, Int_t iact, Double_t step, Double_t *safe) const
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];
452 Double_t r = TMath::Sqrt(r2);
453 // first check if any crossing at all
454 Double_t snxt = TGeoShape::Big();
455 Double_t rdotn = point[0] * dir[0] + point[1] * dir[1] + point[2] * dir[2];
456 Bool_t fullsph = (!TestShapeBit(kGeoThetaSeg) && !TestShapeBit(kGeoPhiSeg)) ? kTRUE : kFALSE;
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) {
465 Bool_t inrmax = kFALSE;
466 Bool_t inrmin = kFALSE;
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
479 if (TGeoShape::IsSameWithinTolerance(fRmin, 0) || rdotn >= 0)
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;
510 Double_t snext = snxt;
511 Double_t st1 = TGeoShape::Big(), st2 = TGeoShape::Big();
512
513 if (TestShapeBit(kGeoThetaSeg)) {
514 if (fTheta1 > 0) {
515 if (TGeoShape::IsSameWithinTolerance(fTheta1, 90)) {
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
523 if (IsPointInside(&ptnew[0], kTRUE, kFALSE, kTRUE))
524 return TMath::Min(snxt, snext);
525 }
526 } else {
527 si = TMath::Sin(fTheta1 * TMath::DegToRad());
528 ci = TMath::Cos(fTheta1 * TMath::DegToRad());
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);
549 Bool_t skip = kFALSE;
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())) {
555 Double_t ddotn =
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;
559 snxt = TGeoShape::Big();
560 }
561 if (snxt < 1E10) {
562 znew = ptnew[2] + snxt * dir[2];
563 if (snxt > 0 && TMath::Abs(znew) < dz) {
564 if (!TestShapeBit(kGeoPhiSeg))
565 st1 = snxt;
566 else {
567 xnew = ptnew[0] + snxt * dir[0];
568 ynew = ptnew[1] + snxt * dir[1];
569 phi0 = TMath::ATan2(ynew, xnew) * TMath::RadToDeg();
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) {
582 if (!TestShapeBit(kGeoPhiSeg))
583 st1 = snxt;
584 else {
585 xnew = ptnew[0] + snxt * dir[0];
586 ynew = ptnew[1] + snxt * dir[1];
587 phi0 = TMath::ATan2(ynew, xnew) * TMath::RadToDeg();
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) {
601 if (TGeoShape::IsSameWithinTolerance(fTheta2, 90)) {
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
609 if (IsPointInside(&ptnew[0], kTRUE, kFALSE, kTRUE))
610 return TMath::Min(snxt, snext);
611 }
612 } else {
613 si = TMath::Sin(fTheta2 * TMath::DegToRad());
614 ci = TMath::Cos(fTheta2 * TMath::DegToRad());
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);
635 Bool_t skip = kFALSE;
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())) {
641 Double_t ddotn =
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;
645 snxt = TGeoShape::Big();
646 }
647 if (snxt < 1E10) {
648 znew = ptnew[2] + snxt * dir[2];
649 if (snxt > 0 && TMath::Abs(znew) < dz) {
650 if (!TestShapeBit(kGeoPhiSeg))
651 st2 = snxt;
652 else {
653 xnew = ptnew[0] + snxt * dir[0];
654 ynew = ptnew[1] + snxt * dir[1];
655 phi0 = TMath::ATan2(ynew, xnew) * TMath::RadToDeg();
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) {
668 if (!TestShapeBit(kGeoPhiSeg))
669 st2 = snxt;
670 else {
671 xnew = ptnew[0] + snxt * dir[0];
672 ynew = ptnew[1] + snxt * dir[1];
673 phi0 = TMath::ATan2(ynew, xnew) * TMath::RadToDeg();
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 }
686 snxt = TMath::Min(st1, st2);
687 snxt = TMath::Min(snxt, snext);
688 // if (snxt<1E20) return snxt;
689 if (TestShapeBit(kGeoPhiSeg)) {
690 Double_t s1 = TMath::Sin(fPhi1 * TMath::DegToRad());
691 Double_t c1 = TMath::Cos(fPhi1 * TMath::DegToRad());
692 Double_t s2 = TMath::Sin(fPhi2 * TMath::DegToRad());
693 Double_t c2 = TMath::Cos(fPhi2 * TMath::DegToRad());
694 Double_t phim = 0.5 * (fPhi1 + fPhi2);
695 Double_t sm = TMath::Sin(phim * TMath::DegToRad());
696 Double_t cm = TMath::Cos(phim * TMath::DegToRad());
697 Double_t s = 0;
698 Double_t safety, un;
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
736Double_t
737TGeoSphere::DistFromInside(const Double_t *point, const Double_t *dir, Int_t iact, Double_t step, Double_t *safe) const
738{
739 Double_t saf[6];
740 Double_t rxy2 = point[0] * point[0] + point[1] * point[1];
741 Double_t rxy = TMath::Sqrt(rxy2);
742 Double_t rad2 = rxy2 + point[2] * point[2];
743 Double_t r = TMath::Sqrt(rad2);
744 Bool_t rzero = kFALSE;
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
755 if (TestShapeBit(kGeoPhiSeg)) {
756 phi = TMath::ATan2(point[1], point[0]) * TMath::RadToDeg();
757 if (phi < 0)
758 phi += 360.;
759 }
760 if (iact < 3 && safe) {
761 saf[0] = (TGeoShape::IsSameWithinTolerance(fRmin, 0)) ? TGeoShape::Big() : r - fRmin;
762 saf[1] = fRmax - r;
763 saf[2] = saf[3] = saf[4] = saf[5] = TGeoShape::Big();
764 if (TestShapeBit(kGeoThetaSeg)) {
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 }
772 if (TestShapeBit(kGeoPhiSeg)) {
773 Double_t dph1 = phi - fPhi1;
774 if (dph1 < 0)
775 dph1 += 360.;
776 if (dph1 <= 90.)
777 saf[4] = rxy * TMath::Sin(dph1 * TMath::DegToRad());
778 Double_t dph2 = fPhi2 - phi;
779 if (dph2 < 0)
780 dph2 += 360.;
781 if (dph2 <= 90.)
782 saf[5] = rxy * TMath::Sin(dph2 * TMath::DegToRad());
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];
798 Double_t sn1 = TGeoShape::Big();
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();
819 if (TestShapeBit(kGeoThetaSeg)) {
820 if (TGeoShape::IsSameWithinTolerance(fTheta1, 90)) {
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];
827 Double_t si = TMath::Sin(fTheta1 * TMath::DegToRad());
828 Double_t ci = TMath::Cos(fTheta1 * TMath::DegToRad());
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())) {
849 Double_t ddotn =
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) {
859 if (!TestShapeBit(kGeoPhiSeg))
860 sn1 = snxt;
861 else {
862 xnew = ptnew[0] + snxt * dir[0];
863 ynew = ptnew[1] + snxt * dir[1];
864 phi0 = TMath::ATan2(ynew, xnew) * TMath::RadToDeg();
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) {
876 if (!TestShapeBit(kGeoPhiSeg))
877 sn1 = snxt;
878 else {
879 xnew = ptnew[0] + snxt * dir[0];
880 ynew = ptnew[1] + snxt * dir[1];
881 phi0 = TMath::ATan2(ynew, xnew) * TMath::RadToDeg();
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 }
894 if (TGeoShape::IsSameWithinTolerance(fTheta2, 90)) {
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];
901 Double_t si = TMath::Sin(fTheta2 * TMath::DegToRad());
902 Double_t ci = TMath::Cos(fTheta2 * TMath::DegToRad());
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())) {
923 Double_t ddotn =
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) {
933 if (!TestShapeBit(kGeoPhiSeg))
934 sn2 = snxt;
935 else {
936 xnew = ptnew[0] + snxt * dir[0];
937 ynew = ptnew[1] + snxt * dir[1];
938 phi0 = TMath::ATan2(ynew, xnew) * TMath::RadToDeg();
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) {
950 if (!TestShapeBit(kGeoPhiSeg))
951 sn2 = snxt;
952 else {
953 xnew = ptnew[0] + snxt * dir[0];
954 ynew = ptnew[1] + snxt * dir[1];
955 phi0 = TMath::ATan2(ynew, xnew) * TMath::RadToDeg();
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 }
969 Double_t st = TMath::Min(sn1, sn2);
970 Double_t sp = TGeoShape::Big();
971 if (TestShapeBit(kGeoPhiSeg)) {
972 Double_t s1 = TMath::Sin(fPhi1 * TMath::DegToRad());
973 Double_t c1 = TMath::Cos(fPhi1 * TMath::DegToRad());
974 Double_t s2 = TMath::Sin(fPhi2 * TMath::DegToRad());
975 Double_t c2 = TMath::Cos(fPhi2 * TMath::DegToRad());
976 Double_t phim = 0.5 * (fPhi1 + fPhi2);
977 Double_t sm = TMath::Sin(phim * TMath::DegToRad());
978 Double_t cm = TMath::Cos(phim * TMath::DegToRad());
979 sp = TGeoShape::DistToPhiMin(point, dir, s1, c1, s2, c2, sm, cm);
980 }
981 Double_t snxt = TMath::Min(sr, st);
982 snxt = TMath::Min(snxt, sp);
983 return snxt;
984}
985
986////////////////////////////////////////////////////////////////////////////////
987/// compute distance to sphere of radius rsph. Direction has to be a unit vector
988
989Double_t TGeoSphere::DistToSphere(const Double_t *point, const Double_t *dir, Double_t rsph, Bool_t check,
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 *
1027TGeoSphere::Divide(TGeoVolume *voldiv, const char *divname, Int_t iaxis, Int_t ndiv, Double_t start, Double_t step)
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);
1039 vmulti = gGeoManager->MakeVolumeMulti(divname, voldiv->GetMedium());
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());
1057 vmulti = gGeoManager->MakeVolumeMulti(divname, 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);
1067 vmulti = gGeoManager->MakeVolumeMulti(divname, voldiv->GetMedium());
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
1086const char *TGeoSphere::GetAxisName(Int_t iaxis) const
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
1099Double_t TGeoSphere::GetAxisRange(Int_t iaxis, Double_t &xlo, Double_t &xhi) const
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
1128void TGeoSphere::GetBoundingCylinder(Double_t *param) const
1129{
1130 Double_t smin = TMath::Sin(fTheta1 * TMath::DegToRad());
1131 Double_t smax = TMath::Sin(fTheta2 * TMath::DegToRad());
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
1156void TGeoSphere::InspectShape() const
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");
1166 TGeoBBox::InspectShape();
1167}
1168
1169////////////////////////////////////////////////////////////////////////////////
1170/// Creates a TBuffer3D describing *this* shape.
1171/// Coordinates are in local reference frame.
1172
1173TBuffer3D *TGeoSphere::MakeBuffer3D() const
1174{
1175 Bool_t full = kTRUE;
1176 if (TestShapeBit(kGeoThetaSeg) || TestShapeBit(kGeoPhiSeg))
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;
1187 if (TestShapeBit(kGeoPhiSeg))
1188 nlong++;
1189
1190 Int_t nbPnts = nlat * nlong + nup + ndown + ncenter;
1191 if (TestShapeBit(kGeoRSeg))
1192 nbPnts *= 2;
1193
1194 Int_t nbSegs = nlat * fNseg + (nlat - 1 + nup + ndown) * nlong; // outer sphere
1195 if (TestShapeBit(kGeoRSeg))
1196 nbSegs *= 2; // inner sphere
1197 if (TestShapeBit(kGeoPhiSeg))
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
1202 if (TestShapeBit(kGeoRSeg))
1203 nbPols *= 2; // inner
1204 if (TestShapeBit(kGeoPhiSeg))
1205 nbPols += 2 * fNz; // 2 phi planes
1206 nbPols += (2 - nup - ndown) * fNseg; // connecting
1207
1208 TBuffer3D *buff =
1209 new TBuffer3D(TBuffer3DTypes::kGeneric, nbPnts, 3 * nbPnts, nbSegs, 3 * nbSegs, nbPols, 6 * nbPols);
1210
1211 if (buff) {
1212 SetPoints(buff->fPnts);
1213 SetSegsAndPols(*buff);
1214 }
1215
1216 return buff;
1217}
1218
1219////////////////////////////////////////////////////////////////////////////////
1220/// Fill TBuffer3D structure for segments and polygons.
1221
1222void TGeoSphere::SetSegsAndPols(TBuffer3D &buff) const
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;
1234 if (TestShapeBit(kGeoPhiSeg))
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)
1267 Int_t indlong = indpar + nlat * fNseg;
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 }
1297 Int_t indparin = inddown + ndown * nlong;
1298 Int_t indlongin = indparin;
1299 Int_t indupin = indparin;
1300 Int_t inddownin = indparin;
1301 Int_t indphi = indparin;
1302 // inner sphere
1303 Int_t indptin = nlat * nlong + nup + ndown;
1304 Int_t iptcenter = indptin;
1305 // nlat*fNseg+(nlat+nup+ndown-1)*nlong
1306 if (TestShapeBit(kGeoRSeg)) {
1307 indlongin = indparin + nlat * fNseg;
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) {
1331 Int_t indupltop = indptin + nlat * nlong;
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) {
1341 Int_t indpdown = indptin + nlat * nlong + nup;
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 }
1348 indphi = inddownin + ndown * nlong;
1349 }
1350 Int_t indtheta = indphi;
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;
1357 if (TestShapeBit(kGeoRSeg))
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;
1365 if (TestShapeBit(kGeoRSeg))
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;
1373 if (TestShapeBit(kGeoRSeg))
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;
1381 if (TestShapeBit(kGeoRSeg))
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;
1392 if (TestShapeBit(kGeoRSeg))
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;
1402 if (TestShapeBit(kGeoRSeg))
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
1605Double_t TGeoSphere::Safety(const Double_t *point, Bool_t in) const
1606{
1607 Double_t r2 = point[0] * point[0] + point[1] * point[1] + point[2] * point[2];
1608 Double_t r = TMath::Sqrt(r2);
1609 Bool_t rzero = kFALSE;
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];
1618 saf[0] = (TGeoShape::IsSameWithinTolerance(fRmin, 0) && !TestShapeBit(kGeoThetaSeg) && !TestShapeBit(kGeoPhiSeg))
1619 ? TGeoShape::Big()
1620 : r - fRmin;
1621 saf[1] = fRmax - r;
1622 saf[2] = saf[3] = TGeoShape::Big();
1623 if (TestShapeBit(kGeoThetaSeg)) {
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 }
1629 Double_t safphi = TGeoShape::Big();
1630 if (TestShapeBit(kGeoPhiSeg))
1631 safphi = TGeoShape::SafetyPhi(point, in, fPhi1, fPhi2);
1632 if (in) {
1633 Double_t safe = saf[TMath::LocMin(4, saf)];
1634 return TMath::Min(safe, safphi);
1635 }
1636 for (Int_t i = 0; i < 4; i++)
1637 saf[i] = -saf[i];
1638 Double_t safe = saf[TMath::LocMax(4, saf)];
1639 if (TestShapeBit(kGeoPhiSeg))
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{
1649 if (TObject::TestBit(kGeoSavePrimitive))
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;
1660 TObject::SetBit(TGeoShape::kGeoSavePrimitive);
1661}
1662
1663////////////////////////////////////////////////////////////////////////////////
1664/// Set spherical segment dimensions.
1665
1666void TGeoSphere::SetSphDimensions(Double_t rmin, Double_t rmax, Double_t theta1, Double_t theta2, Double_t phi1,
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)
1676 SetShapeBit(kGeoRSeg);
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.)
1684 SetShapeBit(kGeoThetaSeg);
1685 fPhi1 = phi1;
1686 if (phi1 < 0)
1687 fPhi1 += 360.;
1688 fPhi2 = phi2;
1689 while (fPhi2 <= fPhi1)
1690 fPhi2 += 360.;
1691 if (!TGeoShape::IsSameWithinTolerance(TMath::Abs(phi2 - phi1), 360))
1692 SetShapeBit(kGeoPhiSeg);
1693}
1694
1695////////////////////////////////////////////////////////////////////////////////
1696/// Set dimensions of the spherical segment starting from a list of parameters.
1697
1698void TGeoSphere::SetDimensions(Double_t *param, Int_t nparam)
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];
1714 SetSphDimensions(rmin, rmax, theta1, theta2, phi1, phi2);
1715}
1716
1717////////////////////////////////////////////////////////////////////////////////
1718/// Set dimensions of the spherical segment starting from a list of parameters.
1719/// Only takes rmin and rmax
1720
1721void TGeoSphere::SetDimensions(Double_t *param)
1722{
1723 SetDimensions(param, 2);
1724}
1725
1726////////////////////////////////////////////////////////////////////////////////
1727/// Set the number of divisions of mesh circles keeping aspect ratio.
1728
1729void TGeoSphere::SetNumberOfDivisions(Int_t p)
1730{
1731 fNseg = p;
1732 Double_t dphi = fPhi2 - fPhi1;
1733 if (dphi < 0)
1734 dphi += 360;
1735 Double_t dtheta = TMath::Abs(fTheta2 - fTheta1);
1736 fNz = Int_t(fNseg * dtheta / dphi) + 1;
1737 if (fNz < 2)
1738 fNz = 2;
1739}
1740
1741////////////////////////////////////////////////////////////////////////////////
1742/// create sphere mesh points
1743
1744void TGeoSphere::SetPoints(Double_t *points) const
1745{
1746 if (!points) {
1747 Error("SetPoints", "Input array is NULL");
1748 return;
1749 }
1750 Bool_t full = kTRUE;
1751 if (TestShapeBit(kGeoThetaSeg) || TestShapeBit(kGeoPhiSeg))
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;
1762 if (TestShapeBit(kGeoPhiSeg))
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;
1768 Double_t phi1 = fPhi1 * TMath::DegToRad();
1769 Double_t phi2 = fPhi2 * TMath::DegToRad();
1770 Double_t dphi = (phi2 - phi1) / fNseg;
1771 Double_t theta1 = fTheta1 * TMath::DegToRad();
1772 Double_t theta2 = fTheta2 * TMath::DegToRad();
1773 Double_t dtheta = (theta2 - theta1) / fNz;
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
1852void TGeoSphere::SetPoints(Float_t *points) const
1853{
1854 if (!points) {
1855 Error("SetPoints", "Input array is NULL");
1856 return;
1857 }
1858 Bool_t full = kTRUE;
1859 if (TestShapeBit(kGeoThetaSeg) || TestShapeBit(kGeoPhiSeg))
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;
1870 if (TestShapeBit(kGeoPhiSeg))
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;
1876 Double_t phi1 = fPhi1 * TMath::DegToRad();
1877 Double_t phi2 = fPhi2 * TMath::DegToRad();
1878 Double_t dphi = (phi2 - phi1) / fNseg;
1879 Double_t theta1 = fTheta1 * TMath::DegToRad();
1880 Double_t theta2 = fTheta2 * TMath::DegToRad();
1881 Double_t dtheta = (theta2 - theta1) / fNz;
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
1960void TGeoSphere::GetMeshNumbers(Int_t &nvert, Int_t &nsegs, Int_t &npols) const
1961{
1962 TGeoSphere *localThis = const_cast<TGeoSphere *>(this);
1963 localThis->SetNumberOfDivisions(gGeoManager->GetNsegments());
1964 Bool_t full = kTRUE;
1965 if (TestShapeBit(kGeoThetaSeg) || TestShapeBit(kGeoPhiSeg))
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;
1976 if (TestShapeBit(kGeoPhiSeg))
1977 nlong++;
1978
1979 nvert = nlat * nlong + nup + ndown + ncenter;
1980 if (TestShapeBit(kGeoRSeg))
1981 nvert *= 2;
1982
1983 nsegs = nlat * fNseg + (nlat - 1 + nup + ndown) * nlong; // outer sphere
1984 if (TestShapeBit(kGeoRSeg))
1985 nsegs *= 2; // inner sphere
1986 if (TestShapeBit(kGeoPhiSeg))
1987 nsegs += 2 * nlat + nup + ndown; // 2 phi planes
1988 nsegs += nlong * (2 - nup - ndown); // connecting cones
1989
1990 npols = fNz * fNseg; // outer
1991 if (TestShapeBit(kGeoRSeg))
1992 npols *= 2; // inner
1993 if (TestShapeBit(kGeoPhiSeg))
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
2001Int_t TGeoSphere::GetNmeshVertices() const
2002{
2003 Bool_t full = kTRUE;
2004 if (TestShapeBit(kGeoThetaSeg) || TestShapeBit(kGeoPhiSeg))
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;
2015 if (TestShapeBit(kGeoPhiSeg))
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;
2021 if (TestShapeBit(kGeoRSeg))
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
2031void TGeoSphere::Sizeof3D() const {}
2032
2033////////////////////////////////////////////////////////////////////////////////
2034/// Fills a static 3D buffer and returns a reference.
2035
2036const TBuffer3D &TGeoSphere::GetBuffer3D(Int_t reqSections, Bool_t localFrame) const
2037{
2038 static TBuffer3DSphere buffer;
2039
2040 TGeoBBox::FillBuffer3D(buffer, reqSections, localFrame);
2041
2042 if (reqSections & TBuffer3D::kShapeSpecific) {
2043 buffer.fRadiusInner = fRmin;
2044 buffer.fRadiusOuter = fRmax;
2045 buffer.fThetaMin = fTheta1;
2046 buffer.fThetaMax = fTheta2;
2047 buffer.fPhiMin = fPhi1;
2048 buffer.fPhiMax = fPhi2;
2049 buffer.SetSectionsValid(TBuffer3D::kShapeSpecific);
2050 }
2051 if (reqSections & TBuffer3D::kRawSizes) {
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;
2057 if (TestShapeBit(kGeoThetaSeg) || TestShapeBit(kGeoPhiSeg))
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;
2068 if (TestShapeBit(kGeoPhiSeg))
2069 nlong++;
2070
2071 Int_t nbPnts = nlat * nlong + nup + ndown + ncenter;
2072 if (TestShapeBit(kGeoRSeg))
2073 nbPnts *= 2;
2074
2075 Int_t nbSegs = nlat * fNseg + (nlat - 1 + nup + ndown) * nlong; // outer sphere
2076 if (TestShapeBit(kGeoRSeg))
2077 nbSegs *= 2; // inner sphere
2078 if (TestShapeBit(kGeoPhiSeg))
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
2083 if (TestShapeBit(kGeoRSeg))
2084 nbPols *= 2; // inner
2085 if (TestShapeBit(kGeoPhiSeg))
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)) {
2090 buffer.SetSectionsValid(TBuffer3D::kRawSizes);
2091 }
2092 }
2093 if ((reqSections & TBuffer3D::kRaw) && buffer.SectionsValid(TBuffer3D::kRawSizes)) {
2094 SetPoints(buffer.fPnts);
2095 if (!buffer.fLocalFrame) {
2096 TransformPoints(buffer.fPnts, buffer.NbPnts());
2097 }
2098 SetSegsAndPols(buffer);
2099 buffer.SetSectionsValid(TBuffer3D::kRaw);
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
2110void TGeoSphere::Contains_v(const Double_t *points, Bool_t *inside, Int_t vecsize) const
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
2121void TGeoSphere::ComputeNormal_v(const Double_t *points, const Double_t *dirs, Double_t *norms, Int_t vecsize)
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
2130void TGeoSphere::DistFromInside_v(const Double_t *points, const Double_t *dirs, Double_t *dists, Int_t vecsize,
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
2140void TGeoSphere::DistFromOutside_v(const Double_t *points, const Double_t *dirs, Double_t *dists, Int_t vecsize,
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
2152void TGeoSphere::Safety_v(const Double_t *points, const Bool_t *inside, Double_t *safe, Int_t vecsize) const
2153{
2154 for (Int_t i = 0; i < vecsize; i++)
2155 safe[i] = Safety(&points[3 * i], inside[i]);
2156}
d
#define d(i)
Definition RSha256.hxx:102
b
#define b(i)
Definition RSha256.hxx:100
c
#define c(i)
Definition RSha256.hxx:101
a
#define a(i)
Definition RSha256.hxx:99
s1
#define s1(x)
Definition RSha256.hxx:91
capacity
std::size_t capacity
Definition RooBatchCompute.cu:105
Int_t
int Int_t
Signed integer 4 bytes (int)
Definition RtypesCore.h:60
Float_t
float Float_t
Float 4 bytes (float)
Definition RtypesCore.h:72
kFALSE
constexpr Bool_t kFALSE
Definition RtypesCore.h:109
Double_t
double Double_t
Double 8 bytes.
Definition RtypesCore.h:74
kTRUE
constexpr Bool_t kTRUE
Definition RtypesCore.h:108
Option_t
const char Option_t
Option string (const char)
Definition RtypesCore.h:81
p
winID h TVirtualViewer3D TVirtualGLPainter p
Definition TGWin32VirtualGLProxy.cxx:52
r
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
Definition TGWin32VirtualXProxy.cxx:171
id
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize id
Definition TGWin32VirtualXProxy.cxx:97
points
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t points
Definition TGWin32VirtualXProxy.cxx:151
name
char name[80]
Definition TGX11.cxx:142
gGeoManager
R__EXTERN TGeoManager * gGeoManager
Definition TGeoManager.h:620
xmin
float xmin
Definition THbookFile.cxx:95
ymin
float ymin
Definition THbookFile.cxx:95
xmax
float xmax
Definition THbookFile.cxx:95
ymax
float ymax
Definition THbookFile.cxx:95
TBuffer3DSphere
Sphere description class - see TBuffer3DTypes for producer classes Supports hollow and cut spheres.
Definition TBuffer3D.h:139
TBuffer3DSphere::fRadiusInner
Double_t fRadiusInner
Definition TBuffer3D.h:153
TBuffer3DSphere::fThetaMin
Double_t fThetaMin
Definition TBuffer3D.h:155
TBuffer3DSphere::fPhiMin
Double_t fPhiMin
Definition TBuffer3D.h:157
TBuffer3DSphere::fThetaMax
Double_t fThetaMax
Definition TBuffer3D.h:156
TBuffer3DSphere::fRadiusOuter
Double_t fRadiusOuter
Definition TBuffer3D.h:154
TBuffer3DSphere::fPhiMax
Double_t fPhiMax
Definition TBuffer3D.h:158
TBuffer3D
Generic 3D primitive description class.
Definition TBuffer3D.h:18
TBuffer3D::fPols
Int_t * fPols
Definition TBuffer3D.h:124
TBuffer3D::NbPnts
UInt_t NbPnts() const
Definition TBuffer3D.h:89
TBuffer3D::SectionsValid
Bool_t SectionsValid(UInt_t mask) const
Definition TBuffer3D.h:76
TBuffer3D::kShapeSpecific
@ kShapeSpecific
Definition TBuffer3D.h:61
TBuffer3D::SetSectionsValid
void SetSectionsValid(UInt_t mask)
Definition TBuffer3D.h:74
TBuffer3D::fSegs
Int_t * fSegs
Definition TBuffer3D.h:123
TBuffer3D::fLocalFrame
Bool_t fLocalFrame
Definition TBuffer3D.h:99
TBuffer3D::SetRawSizes
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
TBuffer3D::fPnts
Double_t * fPnts
Definition TBuffer3D.h:122
TGeoBBox
Box class.
Definition TGeoBBox.h:18
TGeoBBox::FillBuffer3D
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...
Definition TGeoBBox.cxx:1313
TGeoBBox::fDX
Double_t fDX
Definition TGeoBBox.h:21
TGeoBBox::SetBoxDimensions
void SetBoxDimensions(Double_t dx, Double_t dy, Double_t dz, Double_t *origin=nullptr)
Set parameters of the box.
Definition TGeoBBox.cxx:1146
TGeoBBox::DistFromOutside
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
TGeoBBox::fOrigin
Double_t fOrigin[3]
Definition TGeoBBox.h:24
TGeoBBox::InspectShape
void InspectShape() const override
Prints shape parameters.
Definition TGeoBBox.cxx:845
TGeoBBox::fDY
Double_t fDY
Definition TGeoBBox.h:22
TGeoBBox::fDZ
Double_t fDZ
Definition TGeoBBox.h:23
TGeoCone::DistToCone
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
TGeoManager::MakeVolumeMulti
TGeoVolumeMulti * MakeVolumeMulti(const char *name, TGeoMedium *medium)
Make a TGeoVolumeMulti handling a list of volumes.
Definition TGeoManager.cxx:3662
TGeoManager::GetNsegments
Int_t GetNsegments() const
Get number of segments approximating circles.
Definition TGeoManager.cxx:3709
TGeoNodeOffset
Node containing an offset.
Definition TGeoNode.h:185
TGeoPatternFinder
base finder class for patterns. A pattern is specifying a division type
Definition TGeoPatternFinder.h:26
TGeoPatternFinder::SetDivIndex
void SetDivIndex(Int_t index)
Definition TGeoPatternFinder.h:110
TGeoPatternSphPhi
a spherical phi divison pattern
Definition TGeoPatternFinder.h:445
TGeoPatternSphR
a spherical R divison pattern
Definition TGeoPatternFinder.h:395
TGeoPatternSphTheta
a spherical theta divison pattern
Definition TGeoPatternFinder.h:420
TGeoShape
Base abstract class for all shapes.
Definition TGeoShape.h:25
TGeoShape::Big
static Double_t Big()
Definition TGeoShape.h:95
TGeoShape::GetBasicColor
Int_t GetBasicColor() const
Get the basic color (0-7).
Definition TGeoShape.cxx:724
TGeoShape::TransformPoints
void TransformPoints(Double_t *points, UInt_t NbPoints) const
Tranform a set of points (LocalToMaster)
Definition TGeoShape.cxx:596
TGeoShape::SetShapeBit
void SetShapeBit(UInt_t f, Bool_t set)
Equivalent of TObject::SetBit.
Definition TGeoShape.cxx:568
TGeoShape::DistToPhiMin
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.
Definition TGeoShape.cxx:432
TGeoShape::SafetyPhi
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:500
TGeoShape::IsSameWithinTolerance
static Bool_t IsSameWithinTolerance(Double_t a, Double_t b)
Check if two numbers differ with less than a tolerance.
Definition TGeoShape.cxx:344
TGeoShape::GetPointerName
const char * GetPointerName() const
Provide a pointer name containing uid.
Definition TGeoShape.cxx:752
TGeoShape::ShapeDistancetoPrimitive
Int_t ShapeDistancetoPrimitive(Int_t numpoints, Int_t px, Int_t py) const
Returns distance to shape primitive mesh.
Definition TGeoShape.cxx:259
TGeoShape::NormalPhi
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:470
TGeoShape::GetName
const char * GetName() const override
Get the shape name.
Definition TGeoShape.cxx:248
TGeoShape::kGeoSavePrimitive
@ kGeoSavePrimitive
Definition TGeoShape.h:65
TGeoShape::kGeoThetaSeg
@ kGeoThetaSeg
Definition TGeoShape.h:36
TGeoShape::Tolerance
static Double_t Tolerance()
Definition TGeoShape.h:98
TGeoShape::IsCloseToPhi
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,...
Definition TGeoShape.cxx:279
TGeoShape::TestShapeBit
Bool_t TestShapeBit(UInt_t f) const
Definition TGeoShape.h:177
TGeoSphere
TGeoSphere are not just balls having internal and external radii, but sectors of a sphere having defi...
Definition TGeoSphere.h:17
TGeoSphere::DistToSphere
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
Definition TGeoSphere.cxx:989
TGeoSphere::fPhi1
Double_t fPhi1
Definition TGeoSphere.h:25
TGeoSphere::fTheta2
Double_t fTheta2
Definition TGeoSphere.h:24
TGeoSphere::Divide
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,...
Definition TGeoSphere.cxx:1027
TGeoSphere::DistFromInside_v
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...
Definition TGeoSphere.cxx:2130
TGeoSphere::SetPoints
void SetPoints(Double_t *points) const override
create sphere mesh points
Definition TGeoSphere.cxx:1744
TGeoSphere::Contains
Bool_t Contains(const Double_t *point) const override
test if point is inside this sphere check Rmin<=R<=Rmax
Definition TGeoSphere.cxx:386
TGeoSphere::Contains_v
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.
Definition TGeoSphere.cxx:2110
TGeoSphere::SetSphDimensions
void SetSphDimensions(Double_t rmin, Double_t rmax, Double_t theta1, Double_t theta2, Double_t phi1, Double_t phi2)
Set spherical segment dimensions.
Definition TGeoSphere.cxx:1666
TGeoSphere::GetMeshNumbers
void GetMeshNumbers(Int_t &nvert, Int_t &nsegs, Int_t &npols) const override
Returns numbers of vertices, segments and polygons composing the shape mesh.
Definition TGeoSphere.cxx:1960
TGeoSphere::fRmax
Double_t fRmax
Definition TGeoSphere.h:22
TGeoSphere::Safety_v
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.
Definition TGeoSphere.cxx:2152
TGeoSphere::fNseg
Int_t fNseg
Definition TGeoSphere.h:20
TGeoSphere::DistFromOutside_v
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...
Definition TGeoSphere.cxx:2140
TGeoSphere::fRmin
Double_t fRmin
Definition TGeoSphere.h:21
TGeoSphere::TGeoSphere
TGeoSphere()
Default constructor.
Definition TGeoSphere.cxx:66
TGeoSphere::SetDimensions
void SetDimensions(Double_t *param) override
Set dimensions of the spherical segment starting from a list of parameters.
Definition TGeoSphere.cxx:1721
TGeoSphere::MakeBuffer3D
TBuffer3D * MakeBuffer3D() const override
Creates a TBuffer3D describing this shape.
Definition TGeoSphere.cxx:1173
TGeoSphere::InspectShape
void InspectShape() const override
print shape parameters
Definition TGeoSphere.cxx:1156
TGeoSphere::ComputeNormal
void ComputeNormal(const Double_t *point, const Double_t *dir, Double_t *norm) const override
Compute normal to closest surface from POINT.
Definition TGeoSphere.cxx:231
TGeoSphere::SetSegsAndPols
void SetSegsAndPols(TBuffer3D &buff) const override
Fill TBuffer3D structure for segments and polygons.
Definition TGeoSphere.cxx:1222
TGeoSphere::IsPointInside
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.
Definition TGeoSphere.cxx:352
TGeoSphere::Safety
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.
Definition TGeoSphere.cxx:1605
TGeoSphere::ComputeBBox
void ComputeBBox() override
compute bounding box of the sphere
Definition TGeoSphere.cxx:143
TGeoSphere::fTheta1
Double_t fTheta1
Definition TGeoSphere.h:23
TGeoSphere::GetAxisRange
Double_t GetAxisRange(Int_t iaxis, Double_t &xlo, Double_t &xhi) const override
Get range of shape for a given axis.
Definition TGeoSphere.cxx:1099
TGeoSphere::DistFromOutside
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...
Definition TGeoSphere.cxx:435
TGeoSphere::GetBuffer3D
const TBuffer3D & GetBuffer3D(Int_t reqSections, Bool_t localFrame) const override
Fills a static 3D buffer and returns a reference.
Definition TGeoSphere.cxx:2036
TGeoSphere::GetAxisName
const char * GetAxisName(Int_t iaxis) const override
Returns name of axis IAXIS.
Definition TGeoSphere.cxx:1086
TGeoSphere::Capacity
Double_t Capacity() const override
Computes capacity of the shape in [length^3].
Definition TGeoSphere.cxx:129
TGeoSphere::SetNumberOfDivisions
virtual void SetNumberOfDivisions(Int_t p)
Set the number of divisions of mesh circles keeping aspect ratio.
Definition TGeoSphere.cxx:1729
TGeoSphere::ComputeNormal_v
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...
Definition TGeoSphere.cxx:2121
TGeoSphere::GetBoundingCylinder
void GetBoundingCylinder(Double_t *param) const override
Fill vector param[4] with the bounding cylinder parameters.
Definition TGeoSphere.cxx:1128
TGeoSphere::GetNmeshVertices
Int_t GetNmeshVertices() const override
Return number of vertices of the mesh representation.
Definition TGeoSphere.cxx:2001
TGeoSphere::~TGeoSphere
~TGeoSphere() override
destructor
Definition TGeoSphere.cxx:124
TGeoSphere::DistFromInside
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
Definition TGeoSphere.cxx:737
TGeoSphere::IsOnBoundary
Int_t IsOnBoundary(const Double_t *point) const
Check if a point in local sphere coordinates is close to a boundary within shape tolerance.
Definition TGeoSphere.cxx:300
TGeoSphere::SavePrimitive
void SavePrimitive(std::ostream &out, Option_t *option="") override
Save a primitive as a C++ statement(s) on output stream "out".
Definition TGeoSphere.cxx:1647
TGeoSphere::DistancetoPrimitive
Int_t DistancetoPrimitive(Int_t px, Int_t py) override
compute closest distance from point px,py to each corner
Definition TGeoSphere.cxx:422
TGeoSphere::fPhi2
Double_t fPhi2
Definition TGeoSphere.h:26
TGeoSphere::fNz
Int_t fNz
Definition TGeoSphere.h:19
TGeoSphere::Sizeof3D
void Sizeof3D() const override
Definition TGeoSphere.cxx:2031
TGeoVolumeMulti
Volume families.
Definition TGeoVolume.h:269
TGeoVolumeMulti::AddVolume
void AddVolume(TGeoVolume *vol)
Add a volume with valid shape to the list of volumes.
Definition TGeoVolume.cxx:2703
TGeoVolume
TGeoVolume, TGeoVolumeMulti, TGeoVolumeAssembly are the volume classes.
Definition TGeoVolume.h:45
TGeoVolume::AddNodeOffset
void AddNodeOffset(TGeoVolume *vol, Int_t copy_no, Double_t offset=0, Option_t *option="")
Add a division node to the list of nodes.
Definition TGeoVolume.cxx:1134
TGeoVolume::GetMedium
TGeoMedium * GetMedium() const
Definition TGeoVolume.h:178
TGeoVolume::SetFinder
void SetFinder(TGeoPatternFinder *finder)
Definition TGeoVolume.h:247
TGeoVolume::GetNdaughters
Int_t GetNdaughters() const
Definition TGeoVolume.h:371
TGeoVolume::GetNodes
TObjArray * GetNodes()
Definition TGeoVolume.h:172
TObjArray::At
TObject * At(Int_t idx) const override
Definition TObjArray.h:170
TObject::TestBit
R__ALWAYS_INLINE Bool_t TestBit(UInt_t f) const
Definition TObject.h:204
TObject::ClassName
virtual const char * ClassName() const
Returns name of class to which the object belongs.
Definition TObject.cxx:226
TObject::SetBit
void SetBit(UInt_t f, Bool_t set)
Set or unset the user status bits as specified in f.
Definition TObject.cxx:886
TObject::Error
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition TObject.cxx:1096
TString
Basic string class.
Definition TString.h:137
TString::Data
const char * Data() const
Definition TString.h:385
pt
TPaveText * pt
Definition entrylist_figure1.C:7
c1
return c1
Definition legend1.C:41
n
const Int_t n
Definition legend1.C:16
c2
return c2
Definition legend2.C:14
TMath::ACos
Double_t ACos(Double_t)
Returns the principal value of the arc cosine of x, expressed in radians.
Definition TMath.h:645
TMath::LocMin
Long64_t LocMin(Long64_t n, const T *a)
Returns index of array with the minimum element.
Definition TMath.h:995
TMath::Max
Short_t Max(Short_t a, Short_t b)
Returns the largest of a and b.
Definition TMathBase.h:249
TMath::Sign
T1 Sign(T1 a, T2 b)
Returns a value with the magnitude of a and the sign of b.
Definition TMathBase.h:174
TMath::ATan2
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
TMath::LocMax
Long64_t LocMax(Long64_t n, const T *a)
Returns index of array with the maximum element.
Definition TMath.h:1105
TMath::DegToRad
constexpr Double_t DegToRad()
Conversion from degree to radian: .
Definition TMath.h:82
TMath::Sqrt
Double_t Sqrt(Double_t x)
Returns the square root of x.
Definition TMath.h:675
TMath::Min
Short_t Min(Short_t a, Short_t b)
Returns the smallest of a and b.
Definition TMathBase.h:197
TMath::Cos
Double_t Cos(Double_t)
Returns the cosine of an angle of x radians.
Definition TMath.h:607
TMath::Pi
constexpr Double_t Pi()
Definition TMath.h:40
TMath::Sin
Double_t Sin(Double_t)
Returns the sine of an angle of x radians.
Definition TMath.h:601
TMath::RadToDeg
constexpr Double_t RadToDeg()
Conversion from radian to degree: .
Definition TMath.h:75
TMath::Abs
Short_t Abs(Short_t d)
Returns the absolute value of parameter Short_t d.
Definition TMathBase.h:122
th2
auto * th2
Definition textalign.C:18
th1
auto * th1
Definition textalign.C:14