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