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