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