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