Logo ROOT  

Reference Guide
 
Loading...
Searching...
No Matches
TGeoSphere.cxx
Go to the documentation of this file.
1// @(#)root/geom:$Id$
2// Author: Andrei Gheata 31/01/02
3// TGeoSphere::Contains() DistFromOutside/Out() implemented by Mihaela Gheata
4
5/*************************************************************************
6 * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers. *
7 * All rights reserved. *
8 * *
9 * For the licensing terms see $ROOTSYS/LICENSE. *
10 * For the list of contributors see $ROOTSYS/README/CREDITS. *
11 *************************************************************************/
12
13/** \class TGeoSphere
14\ingroup Shapes_classes
15
16TGeoSphere are not just balls having internal and external
17radii, but sectors of a sphere having defined theta and phi ranges. The
18TGeoSphere class has the following constructor.
19
20~~~{.cpp}
21TGeoSphere(Double_t rmin,Double_t rmax,Double_t theta1,
22Double_t theta2,Double_t phi1, Double_t phi2);
23~~~
24
25Begin_Macro
26{
27 TCanvas *c = new TCanvas("c", "c",0,0,600,600);
28 new TGeoManager("sphere", "poza7");
29 TGeoMaterial *mat = new TGeoMaterial("Al", 26.98,13,2.7);
30 TGeoMedium *med = new TGeoMedium("MED",1,mat);
31 TGeoVolume *top = gGeoManager->MakeBox("TOP",med,100,100,100);
32 gGeoManager->SetTopVolume(top);
33 TGeoVolume *vol = gGeoManager->MakeSphere("SPHERE",med, 30,40,60,120,30,240);
34 vol->SetLineWidth(2);
35 top->AddNode(vol,1);
36 gGeoManager->CloseGeometry();
37 gGeoManager->SetNsegments(30);
38 top->Draw();
39 TView *view = gPad->GetView();
40 view->ShowAxis();
41}
42End_Macro
43
44 - `rmin: ` internal radius of the spherical sector
45 - `rmax:` external radius
46 - `theta1:` starting theta value [0, 180) in degrees
47 - `theta2:` ending theta value (0, 180] in degrees (`theta1<theta2`)
48
49*/
50
51#include <iostream>
52
53#include "TGeoCone.h"
54#include "TGeoManager.h"
55#include "TGeoVolume.h"
56#include "TVirtualGeoPainter.h"
57#include "TGeoSphere.h"
58#include "TBuffer3D.h"
59#include "TBuffer3DTypes.h"
60#include "TMath.h"
61
62ClassImp(TGeoSphere);
63
64////////////////////////////////////////////////////////////////////////////////
65/// Default constructor
66
67TGeoSphere::TGeoSphere()
68{
69 SetShapeBit(TGeoShape::kGeoSph);
70 fNz = 0;
71 fNseg = 0;
72 fRmin = 0.0;
73 fRmax = 0.0;
74 fTheta1 = 0.0;
75 fTheta2 = 180.0;
76 fPhi1 = 0.0;
77 fPhi2 = 360.0;
78}
79
80////////////////////////////////////////////////////////////////////////////////
81/// Default constructor specifying minimum and maximum radius
82
83TGeoSphere::TGeoSphere(Double_t rmin, Double_t rmax, Double_t theta1,
84 Double_t theta2, Double_t phi1, Double_t phi2)
85 :TGeoBBox(0, 0, 0)
86{
87 SetShapeBit(TGeoShape::kGeoSph);
88 SetSphDimensions(rmin, rmax, theta1, theta2, phi1, phi2);
89 ComputeBBox();
90 SetNumberOfDivisions(20);
91}
92
93////////////////////////////////////////////////////////////////////////////////
94/// Default constructor specifying minimum and maximum radius
95
96TGeoSphere::TGeoSphere(const char *name, Double_t rmin, Double_t rmax, Double_t theta1,
97 Double_t theta2, Double_t phi1, Double_t phi2)
98 :TGeoBBox(name, 0, 0, 0)
99{
100 SetShapeBit(TGeoShape::kGeoSph);
101 SetSphDimensions(rmin, rmax, theta1, theta2, phi1, phi2);
102 ComputeBBox();
103 SetNumberOfDivisions(20);
104}
105
106////////////////////////////////////////////////////////////////////////////////
107/// Default constructor specifying minimum and maximum radius
108/// param[0] = Rmin
109/// param[1] = Rmax
110/// param[2] = theta1
111/// param[3] = theta2
112/// param[4] = phi1
113/// param[5] = phi2
114
115TGeoSphere::TGeoSphere(Double_t *param, Int_t nparam)
116 :TGeoBBox(0, 0, 0)
117{
118 SetShapeBit(TGeoShape::kGeoSph);
119 SetDimensions(param, nparam);
120 ComputeBBox();
121 SetNumberOfDivisions(20);
122}
123
124////////////////////////////////////////////////////////////////////////////////
125/// destructor
126
127TGeoSphere::~TGeoSphere()
128{
129}
130
131////////////////////////////////////////////////////////////////////////////////
132/// Computes capacity of the shape in [length^3]
133
134Double_t TGeoSphere::Capacity() const
135{
136 Double_t th1 = fTheta1*TMath::DegToRad();
137 Double_t th2 = fTheta2*TMath::DegToRad();
138 Double_t ph1 = fPhi1*TMath::DegToRad();
139 Double_t ph2 = fPhi2*TMath::DegToRad();
140 Double_t capacity = (1./3.)*(fRmax*fRmax*fRmax-fRmin*fRmin*fRmin)*
141 TMath::Abs(TMath::Cos(th1)-TMath::Cos(th2))*
142 TMath::Abs(ph2-ph1);
143 return capacity;
144}
145
146////////////////////////////////////////////////////////////////////////////////
147/// compute bounding box of the sphere
148
149void TGeoSphere::ComputeBBox()
150{
151 if (TGeoShape::IsSameWithinTolerance(TMath::Abs(fTheta2-fTheta1),180)) {
152 if (TGeoShape::IsSameWithinTolerance(TMath::Abs(fPhi2-fPhi1),360)) {
153 TGeoBBox::SetBoxDimensions(fRmax, fRmax, fRmax);
154 memset(fOrigin, 0, 3*sizeof(Double_t));
155 return;
156 }
157 }
158 Double_t st1 = TMath::Sin(fTheta1*TMath::DegToRad());
159 Double_t st2 = TMath::Sin(fTheta2*TMath::DegToRad());
160 Double_t r1min, r1max, r2min, r2max, rmin, rmax;
161 r1min = TMath::Min(fRmax*st1, fRmax*st2);
162 r1max = TMath::Max(fRmax*st1, fRmax*st2);
163 r2min = TMath::Min(fRmin*st1, fRmin*st2);
164 r2max = TMath::Max(fRmin*st1, fRmin*st2);
165 if (((fTheta1<=90) && (fTheta2>=90)) || ((fTheta2<=90) && (fTheta1>=90))) {
166 r1max = fRmax;
167 r2max = fRmin;
168 }
169 rmin = TMath::Min(r1min, r2min);
170 rmax = TMath::Max(r1max, r2max);
171
172 Double_t xc[4];
173 Double_t yc[4];
174 xc[0] = rmax*TMath::Cos(fPhi1*TMath::DegToRad());
175 yc[0] = rmax*TMath::Sin(fPhi1*TMath::DegToRad());
176 xc[1] = rmax*TMath::Cos(fPhi2*TMath::DegToRad());
177 yc[1] = rmax*TMath::Sin(fPhi2*TMath::DegToRad());
178 xc[2] = rmin*TMath::Cos(fPhi1*TMath::DegToRad());
179 yc[2] = rmin*TMath::Sin(fPhi1*TMath::DegToRad());
180 xc[3] = rmin*TMath::Cos(fPhi2*TMath::DegToRad());
181 yc[3] = rmin*TMath::Sin(fPhi2*TMath::DegToRad());
182
183 Double_t xmin = xc[TMath::LocMin(4, &xc[0])];
184 Double_t xmax = xc[TMath::LocMax(4, &xc[0])];
185 Double_t ymin = yc[TMath::LocMin(4, &yc[0])];
186 Double_t ymax = yc[TMath::LocMax(4, &yc[0])];
187 Double_t dp = fPhi2-fPhi1;
188 if (dp<0) dp+=360;
189 Double_t ddp = -fPhi1;
190 if (ddp<0) ddp+= 360;
191 if (ddp>360) ddp-=360;
192 if (ddp<=dp) xmax = rmax;
193 ddp = 90-fPhi1;
194 if (ddp<0) ddp+= 360;
195 if (ddp>360) ddp-=360;
196 if (ddp<=dp) ymax = rmax;
197 ddp = 180-fPhi1;
198 if (ddp<0) ddp+= 360;
199 if (ddp>360) ddp-=360;
200 if (ddp<=dp) xmin = -rmax;
201 ddp = 270-fPhi1;
202 if (ddp<0) ddp+= 360;
203 if (ddp>360) ddp-=360;
204 if (ddp<=dp) ymin = -rmax;
205 xc[0] = fRmax*TMath::Cos(fTheta1*TMath::DegToRad());
206 xc[1] = fRmax*TMath::Cos(fTheta2*TMath::DegToRad());
207 xc[2] = fRmin*TMath::Cos(fTheta1*TMath::DegToRad());
208 xc[3] = fRmin*TMath::Cos(fTheta2*TMath::DegToRad());
209 Double_t zmin = xc[TMath::LocMin(4, &xc[0])];
210 Double_t zmax = xc[TMath::LocMax(4, &xc[0])];
211
212
213 fOrigin[0] = (xmax+xmin)/2;
214 fOrigin[1] = (ymax+ymin)/2;
215 fOrigin[2] = (zmax+zmin)/2;;
216 fDX = (xmax-xmin)/2;
217 fDY = (ymax-ymin)/2;
218 fDZ = (zmax-zmin)/2;
219}
220
221////////////////////////////////////////////////////////////////////////////////
222/// Compute normal to closest surface from POINT.
223
224void TGeoSphere::ComputeNormal(const Double_t *point, const Double_t *dir, Double_t *norm)
225{
226 Double_t rxy2 = point[0]*point[0]+point[1]*point[1];
227 Double_t r2 = rxy2+point[2]*point[2];
228 Double_t r=TMath::Sqrt(r2);
229 Bool_t rzero=kFALSE;
230 if (r<=1E-20) rzero=kTRUE;
231 //localize theta
232 Double_t phi=0;
233 Double_t th=0.;
234 if (!rzero) th = TMath::ACos(point[2]/r);
235
236 //localize phi
237 phi=TMath::ATan2(point[1], point[0]);
238
239 Double_t saf[4];
240 saf[0]=(TGeoShape::IsSameWithinTolerance(fRmin,0) && !TestShapeBit(kGeoThetaSeg) && !TestShapeBit(kGeoPhiSeg))?TGeoShape::Big():TMath::Abs(r-fRmin);
241 saf[1]=TMath::Abs(fRmax-r);
242 saf[2]=saf[3]= TGeoShape::Big();
243 if (TestShapeBit(kGeoThetaSeg)) {
244 if (fTheta1>0) {
245 saf[2] = r*TMath::Abs(TMath::Sin(th-fTheta1*TMath::DegToRad()));
246 }
247 if (fTheta2<180) {
248 saf[3] = r*TMath::Abs(TMath::Sin(fTheta2*TMath::DegToRad()-th));
249 }
250 }
251 Int_t i = TMath::LocMin(4,saf);
252 if (TestShapeBit(kGeoPhiSeg)) {
253 Double_t c1 = TMath::Cos(fPhi1*TMath::DegToRad());
254 Double_t s1 = TMath::Sin(fPhi1*TMath::DegToRad());
255 Double_t c2 = TMath::Cos(fPhi2*TMath::DegToRad());
256 Double_t s2 = TMath::Sin(fPhi2*TMath::DegToRad());
257 if (TGeoShape::IsCloseToPhi(saf[i], point,c1,s1,c2,s2)) {
258 TGeoShape::NormalPhi(point,dir,norm,c1,s1,c2,s2);
259 return;
260 }
261 }
262 if (i>1) {
263 if (i==2) th=(fTheta1<90)?(fTheta1+90):(fTheta1-90);
264 else th=(fTheta2<90)?(fTheta2+90):(fTheta2-90);
265 th *= TMath::DegToRad();
266 }
267
268 norm[0] = TMath::Sin(th)*TMath::Cos(phi);
269 norm[1] = TMath::Sin(th)*TMath::Sin(phi);
270 norm[2] = TMath::Cos(th);
271 if (norm[0]*dir[0]+norm[1]*dir[1]+norm[2]*dir[2]<0) {
272 norm[0] = -norm[0];
273 norm[1] = -norm[1];
274 norm[2] = -norm[2];
275 }
276}
277
278////////////////////////////////////////////////////////////////////////////////
279/// Check if a point in local sphere coordinates is close to a boundary within
280/// shape tolerance. Return values:
281/// - 0 - not close to boundary
282/// - 1 - close to Rmin boundary
283/// - 2 - close to Rmax boundary
284/// - 3,4 - close to phi1/phi2 boundary
285/// - 5,6 - close to theta1/theta2 boundary
286
287Int_t TGeoSphere::IsOnBoundary(const Double_t *point) const
288{
289 Int_t icode = 0;
290 Double_t tol = TGeoShape::Tolerance();
291 Double_t r2 = point[0]*point[0]+point[1]*point[1]+point[2]*point[2];
292 Double_t drsqout = r2-fRmax*fRmax;
293 // Test if point is on fRmax boundary
294 if (TMath::Abs(drsqout)<2.*fRmax*tol) return 2;
295 Double_t drsqin = r2;
296 // Test if point is on fRmin boundary
297 if (TestShapeBit(kGeoRSeg)) {
298 drsqin -= fRmin*fRmin;
299 if (TMath::Abs(drsqin)<2.*fRmin*tol) return 1;
300 }
301 if (TestShapeBit(kGeoPhiSeg)) {
302 Double_t phi = TMath::ATan2(point[1], point[0]);
303 if (phi<0) phi+=2*TMath::Pi();
304 Double_t phi1 = fPhi1*TMath::DegToRad();
305 Double_t phi2 = fPhi2*TMath::DegToRad();
306 Double_t ddp = phi-phi1;
307 if (r2*ddp*ddp < tol*tol) return 3;
308 ddp = phi - phi2;
309 if (r2*ddp*ddp < tol*tol) return 4;
310 }
311 if (TestShapeBit(kGeoThetaSeg)) {
312 Double_t r = TMath::Sqrt(r2);
313 Double_t theta = TMath::ACos(point[2]/r2);
314 Double_t theta1 = fTheta1*TMath::DegToRad();
315 Double_t theta2 = fTheta2*TMath::DegToRad();
316 Double_t ddt;
317 if (fTheta1>0) {
318 ddt = TMath::Abs(theta-theta1);
319 if (r*ddt < tol) return 5;
320 }
321 if (fTheta2<180) {
322 ddt = TMath::Abs(theta-theta2);
323 if (r*ddt < tol) return 6;
324 }
325 }
326 return icode;
327}
328
329////////////////////////////////////////////////////////////////////////////////
330/// Check if a point is inside radius/theta/phi ranges for the spherical sector.
331
332Bool_t TGeoSphere::IsPointInside(const Double_t *point, Bool_t checkR, Bool_t checkTh, Bool_t checkPh) const
333{
334 Double_t r2 = point[0]*point[0]+point[1]*point[1]+point[2]*point[2];
335 if (checkR) {
336 if (TestShapeBit(kGeoRSeg) && (r2<fRmin*fRmin)) return kFALSE;
337 if (r2>fRmax*fRmax) return kFALSE;
338 }
339 if (r2<1E-20) return kTRUE;
340 if (checkPh && TestShapeBit(kGeoPhiSeg)) {
341 Double_t phi = TMath::ATan2(point[1], point[0]) * TMath::RadToDeg();
342 while (phi < fPhi1) phi+=360.;
343 Double_t dphi = fPhi2 -fPhi1;
344 Double_t ddp = phi - fPhi1;
345 if (ddp > dphi) return kFALSE;
346 }
347 if (checkTh && TestShapeBit(kGeoThetaSeg)) {
348 r2=TMath::Sqrt(r2);
349 // check theta range
350 Double_t theta = TMath::ACos(point[2]/r2)*TMath::RadToDeg();
351 if ((theta<fTheta1) || (theta>fTheta2)) return kFALSE;
352 }
353 return kTRUE;
354}
355
356////////////////////////////////////////////////////////////////////////////////
357/// test if point is inside this sphere
358/// check Rmin<=R<=Rmax
359
360Bool_t TGeoSphere::Contains(const Double_t *point) const
361{
362 Double_t r2=point[0]*point[0]+point[1]*point[1]+point[2]*point[2];
363 if (TestShapeBit(kGeoRSeg) && (r2<fRmin*fRmin)) return kFALSE;
364 if (r2>fRmax*fRmax) return kFALSE;
365 if (r2<1E-20) return kTRUE;
366 // check phi range
367 if (TestShapeBit(kGeoPhiSeg)) {
368 Double_t phi = TMath::ATan2(point[1], point[0]) * TMath::RadToDeg();
369 if (phi < 0 ) phi+=360.;
370 Double_t dphi = fPhi2 -fPhi1;
371 if (dphi < 0) dphi+=360.;
372 Double_t ddp = phi - fPhi1;
373 if (ddp < 0) ddp += 360.;
374 if (ddp > dphi) return kFALSE;
375 }
376 if (TestShapeBit(kGeoThetaSeg)) {
377 r2=TMath::Sqrt(r2);
378 // check theta range
379 Double_t theta = TMath::ACos(point[2]/r2)*TMath::RadToDeg();
380 if ((theta<fTheta1) || (theta>fTheta2)) return kFALSE;
381 }
382 return kTRUE;
383}
384
385////////////////////////////////////////////////////////////////////////////////
386/// compute closest distance from point px,py to each corner
387
388Int_t TGeoSphere::DistancetoPrimitive(Int_t px, Int_t py)
389{
390 Int_t n = fNseg+1;
391 Int_t nz = fNz+1;
392 const Int_t numPoints = 2*n*nz;
393 return ShapeDistancetoPrimitive(numPoints, px, py);
394}
395
396////////////////////////////////////////////////////////////////////////////////
397/// compute distance from outside point to surface of the sphere
398/// Check if the bounding box is crossed within the requested distance
399
400Double_t TGeoSphere::DistFromOutside(const Double_t *point, const Double_t *dir, Int_t iact, Double_t step, Double_t *safe) const
401{
402 Double_t sdist = TGeoBBox::DistFromOutside(point,dir, fDX, fDY, fDZ, fOrigin, step);
403 if (sdist>=step) return TGeoShape::Big();
404 Double_t saf[6];
405 Double_t r1,r2,z1,z2,dz,si,ci;
406 Double_t rxy2 = point[0]*point[0]+point[1]*point[1];
407 Double_t rxy = TMath::Sqrt(rxy2);
408 r2 = rxy2+point[2]*point[2];
409 Double_t r=TMath::Sqrt(r2);
410 Bool_t rzero=kFALSE;
411 Double_t phi=0;
412 if (r<1E-20) rzero=kTRUE;
413 //localize theta
414 Double_t th=0.;
415 if (TestShapeBit(kGeoThetaSeg) && (!rzero)) {
416 th = TMath::ACos(point[2]/r)*TMath::RadToDeg();
417 }
418 //localize phi
419 if (TestShapeBit(kGeoPhiSeg)) {
420 phi=TMath::ATan2(point[1], point[0])*TMath::RadToDeg();
421 if (phi<0) phi+=360.;
422 }
423 if (iact<3 && safe) {
424 saf[0]=(r<fRmin)?fRmin-r:TGeoShape::Big();
425 saf[1]=(r>fRmax)?(r-fRmax):TGeoShape::Big();
426 saf[2]=saf[3]=saf[4]=saf[5]= TGeoShape::Big();
427 if (TestShapeBit(kGeoThetaSeg)) {
428 if (th < fTheta1) {
429 saf[2] = r*TMath::Sin((fTheta1-th)*TMath::DegToRad());
430 }
431 if (th > fTheta2) {
432 saf[3] = r*TMath::Sin((th-fTheta2)*TMath::DegToRad());
433 }
434 }
435 if (TestShapeBit(kGeoPhiSeg)) {
436 Double_t dph1=phi-fPhi1;
437 if (dph1<0) dph1+=360.;
438 if (dph1<=90.) saf[4]=rxy*TMath::Sin(dph1*TMath::DegToRad());
439 Double_t dph2=fPhi2-phi;
440 if (dph2<0) dph2+=360.;
441 if (dph2>90.) saf[5]=rxy*TMath::Sin(dph2*TMath::DegToRad());
442 }
443 *safe = saf[TMath::LocMin(6, &saf[0])];
444 if (iact==0) return TGeoShape::Big();
445 if (iact==1 && step<*safe) return TGeoShape::Big();
446 }
447 // compute distance to shape
448 // first check if any crossing at all
449 Double_t snxt = TGeoShape::Big();
450 Double_t rdotn = point[0]*dir[0]+point[1]*dir[1]+point[2]*dir[2];
451 Bool_t fullsph = (!TestShapeBit(kGeoThetaSeg) && !TestShapeBit(kGeoPhiSeg))?kTRUE:kFALSE;
452 if (r>fRmax) {
453 Double_t b = rdotn;
454 Double_t c = r2-fRmax*fRmax;
455 Double_t d=b*b-c;
456 if (d<0) return TGeoShape::Big();
457 }
458 if (fullsph) {
459 Bool_t inrmax = kFALSE;
460 Bool_t inrmin = kFALSE;
461 if (r<=fRmax+TGeoShape::Tolerance()) inrmax = kTRUE;
462 if (r>=fRmin-TGeoShape::Tolerance()) inrmin = kTRUE;
463 if (inrmax && inrmin) {
464 if ((fRmax-r) < (r-fRmin)) {
465 // close to Rmax
466 if (rdotn>=0) return TGeoShape::Big();
467 return 0.0; // already in
468 }
469 // close to Rmin
470 if (TGeoShape::IsSameWithinTolerance(fRmin,0) || rdotn>=0) return 0.0;
471 // check second crossing of Rmin
472 return DistToSphere(point, dir, fRmin, kFALSE, kFALSE);
473 }
474 }
475
476 // do rmin, rmax, checking phi and theta ranges
477 if (r<fRmin) {
478 // check first cross of rmin
479 snxt = DistToSphere(point, dir, fRmin, kTRUE);
480 if (snxt<1E20) return snxt;
481 } else {
482 if (r>fRmax) {
483 // point outside rmax, check first cross of rmax
484 snxt = DistToSphere(point, dir, fRmax, kTRUE);
485 if (snxt<1E20) return snxt;
486 // now check second crossing of rmin
487 if (fRmin>0) snxt = DistToSphere(point, dir, fRmin, kTRUE, kFALSE);
488 } else {
489 // point between rmin and rmax, check second cross of rmin
490 if (fRmin>0) snxt = DistToSphere(point, dir, fRmin, kTRUE, kFALSE);
491 }
492 }
493 // check theta conical surfaces
494 Double_t ptnew[3];
495 Double_t b,delta,xnew,ynew,znew, phi0, ddp;
496 Double_t snext = snxt;
497 Double_t st1=TGeoShape::Big(), st2=TGeoShape::Big();
498
499 if (TestShapeBit(kGeoThetaSeg)) {
500 if (fTheta1>0) {
501 if (TGeoShape::IsSameWithinTolerance(fTheta1,90)) {
502 // surface is a plane
503 if (point[2]*dir[2]<0) {
504 snxt = -point[2]/dir[2];
505 ptnew[0] = point[0]+snxt*dir[0];
506 ptnew[1] = point[1]+snxt*dir[1];
507 ptnew[2] = 0;
508 // check range
509 if (IsPointInside(&ptnew[0], kTRUE, kFALSE, kTRUE)) return TMath::Min(snxt,snext);
510 }
511 } else {
512 si = TMath::Sin(fTheta1*TMath::DegToRad());
513 ci = TMath::Cos(fTheta1*TMath::DegToRad());
514 if (ci>0) {
515 r1 = fRmin*si;
516 z1 = fRmin*ci;
517 r2 = fRmax*si;
518 z2 = fRmax*ci;
519 } else {
520 r1 = fRmax*si;
521 z1 = fRmax*ci;
522 r2 = fRmin*si;
523 z2 = fRmin*ci;
524 }
525 dz = 0.5*(z2-z1);
526 ptnew[0] = point[0];
527 ptnew[1] = point[1];
528 ptnew[2] = point[2]-0.5*(z1+z2);
529 Double_t zinv = 1./dz;
530 Double_t rin = 0.5*(r1+r2+(r2-r1)*ptnew[2]*zinv);
531 // Protection in case point is outside
532 Double_t sigz = TMath::Sign(1.,point[2]);
533 TGeoCone::DistToCone(ptnew, dir, dz, r1, r2, b, delta);
534 Bool_t skip = kFALSE;
535 if (delta<0) skip = kTRUE;
536 if (!skip) {
537 snxt = -b-delta;
538 if (sigz*ci>0 && sigz*rxy2 > sigz*rin*(rin-sigz*TGeoShape::Tolerance())) {
539 Double_t ddotn = ptnew[0]*dir[0]+ptnew[1]*dir[1]+0.5*(r1-r2)*dir[2]*zinv*TMath::Sqrt(rxy2);
540 if (sigz*ddotn>=0 || -b+delta<1.E-9) skip = kTRUE;
541 snxt = TGeoShape::Big();
542 }
543 if (snxt<1E10) {
544 znew = ptnew[2]+snxt*dir[2];
545 if (snxt>0 && TMath::Abs(znew)<dz) {
546 if (!TestShapeBit(kGeoPhiSeg)) st1=snxt;
547 else {
548 xnew = ptnew[0]+snxt*dir[0];
549 ynew = ptnew[1]+snxt*dir[1];
550 phi0 = TMath::ATan2(ynew, xnew) * TMath::RadToDeg();
551 ddp = phi0-fPhi1;
552 while (ddp<0) ddp+=360.;
553 if (ddp<=fPhi2-fPhi1) st1 = snxt;
554 }
555 }
556 }
557 if (!skip && st1>1E10) {
558 snxt = -b+delta;
559 znew = ptnew[2]+snxt*dir[2];
560 if (snxt>0 && TMath::Abs(znew)<dz) {
561 if (!TestShapeBit(kGeoPhiSeg)) st1=snxt;
562 else {
563 xnew = ptnew[0]+snxt*dir[0];
564 ynew = ptnew[1]+snxt*dir[1];
565 phi0 = TMath::ATan2(ynew, xnew) * TMath::RadToDeg();
566 ddp = phi0-fPhi1;
567 while (ddp<0) ddp+=360.;
568 if (ddp<=fPhi2-fPhi1) st1 = snxt;
569 }
570 }
571 }
572 }
573 }
574 }
575
576 if (fTheta2<180) {
577 if (TGeoShape::IsSameWithinTolerance(fTheta2,90)) {
578 // surface is a plane
579 if (point[2]*dir[2]<0) {
580 snxt = -point[2]/dir[2];
581 ptnew[0] = point[0]+snxt*dir[0];
582 ptnew[1] = point[1]+snxt*dir[1];
583 ptnew[2] = 0;
584 // check range
585 if (IsPointInside(&ptnew[0], kTRUE, kFALSE, kTRUE)) return TMath::Min(snxt,snext);
586 }
587 } else {
588 si = TMath::Sin(fTheta2*TMath::DegToRad());
589 ci = TMath::Cos(fTheta2*TMath::DegToRad());
590 if (ci>0) {
591 r1 = fRmin*si;
592 z1 = fRmin*ci;
593 r2 = fRmax*si;
594 z2 = fRmax*ci;
595 } else {
596 r1 = fRmax*si;
597 z1 = fRmax*ci;
598 r2 = fRmin*si;
599 z2 = fRmin*ci;
600 }
601 dz = 0.5*(z2-z1);
602 ptnew[0] = point[0];
603 ptnew[1] = point[1];
604 ptnew[2] = point[2]-0.5*(z1+z2);
605 Double_t zinv = 1./dz;
606 Double_t rin = 0.5*(r1+r2+(r2-r1)*ptnew[2]*zinv);
607 // Protection in case point is outside
608 Double_t sigz = TMath::Sign(1.,point[2]);
609 TGeoCone::DistToCone(ptnew, dir, dz, r1, r2, b, delta);
610 Bool_t skip = kFALSE;
611 if (delta<0) skip = kTRUE;
612 if (!skip) {
613 snxt = -b-delta;
614 if (sigz*ci>0 && sigz*rxy2 < sigz*rin*(rin+sigz*TGeoShape::Tolerance())) {
615 Double_t ddotn = ptnew[0]*dir[0]+ptnew[1]*dir[1]+0.5*(r1-r2)*dir[2]*zinv*TMath::Sqrt(rxy2);
616 if (sigz*ddotn<=0 || -b+delta<1.E-9) skip = kTRUE;
617 snxt = TGeoShape::Big();
618 }
619 if (snxt<1E10) {
620 znew = ptnew[2]+snxt*dir[2];
621 if (snxt>0 && TMath::Abs(znew)<dz) {
622 if (!TestShapeBit(kGeoPhiSeg)) st2=snxt;
623 else {
624 xnew = ptnew[0]+snxt*dir[0];
625 ynew = ptnew[1]+snxt*dir[1];
626 phi0 = TMath::ATan2(ynew, xnew) * TMath::RadToDeg();
627 ddp = phi0-fPhi1;
628 while (ddp<0) ddp+=360.;
629 if (ddp<=fPhi2-fPhi1) st2 = snxt;
630 }
631 }
632 }
633 if (!skip && st2>1E10) {
634 snxt = -b+delta;
635 znew = ptnew[2]+snxt*dir[2];
636 if (snxt>0 && TMath::Abs(znew)<dz) {
637 if (!TestShapeBit(kGeoPhiSeg)) st2=snxt;
638 else {
639 xnew = ptnew[0]+snxt*dir[0];
640 ynew = ptnew[1]+snxt*dir[1];
641 phi0 = TMath::ATan2(ynew, xnew) * TMath::RadToDeg();
642 ddp = phi0-fPhi1;
643 while (ddp<0) ddp+=360.;
644 if (ddp<=fPhi2-fPhi1) st2 = snxt;
645 }
646 }
647 }
648 }
649 }
650 }
651 }
652 snxt = TMath::Min(st1, st2);
653 snxt = TMath::Min(snxt,snext);
654// if (snxt<1E20) return snxt;
655 if (TestShapeBit(kGeoPhiSeg)) {
656 Double_t s1 = TMath::Sin(fPhi1*TMath::DegToRad());
657 Double_t c1 = TMath::Cos(fPhi1*TMath::DegToRad());
658 Double_t s2 = TMath::Sin(fPhi2*TMath::DegToRad());
659 Double_t c2 = TMath::Cos(fPhi2*TMath::DegToRad());
660 Double_t phim = 0.5*(fPhi1+fPhi2);
661 Double_t sm = TMath::Sin(phim*TMath::DegToRad());
662 Double_t cm = TMath::Cos(phim*TMath::DegToRad());
663 Double_t s=0;
664 Double_t safety, un;
665 safety = point[0]*s1-point[1]*c1;
666 if (safety>0) {
667 un = dir[0]*s1-dir[1]*c1;
668 if (un<0) {
669 s=-safety/un;
670 ptnew[0] = point[0]+s*dir[0];
671 ptnew[1] = point[1]+s*dir[1];
672 ptnew[2] = point[2]+s*dir[2];
673 if ((ptnew[1]*cm-ptnew[0]*sm)<=0) {
674 Double_t sfi1 = s;
675 if (IsPointInside(&ptnew[0], kTRUE, kTRUE, kFALSE) && sfi1<snxt) return sfi1;
676 }
677 }
678 }
679 safety = -point[0]*s2+point[1]*c2;
680 if (safety>0) {
681 un = -dir[0]*s2+dir[1]*c2;
682 if (un<0) {
683 s=-safety/un;
684 ptnew[0] = point[0]+s*dir[0];
685 ptnew[1] = point[1]+s*dir[1];
686 ptnew[2] = point[2]+s*dir[2];
687 if ((ptnew[1]*cm-ptnew[0]*sm)>=0) {
688 Double_t sfi2 = s;
689 if (IsPointInside(&ptnew[0], kTRUE, kTRUE, kFALSE) && sfi2<snxt) return sfi2;
690 }
691 }
692 }
693 }
694 return snxt;
695}
696
697////////////////////////////////////////////////////////////////////////////////
698/// compute distance from inside point to surface of the sphere
699
700Double_t TGeoSphere::DistFromInside(const Double_t *point, const Double_t *dir, Int_t iact, Double_t step, Double_t *safe) const
701{
702 Double_t saf[6];
703 Double_t rxy2 = point[0]*point[0]+point[1]*point[1];
704 Double_t rxy = TMath::Sqrt(rxy2);
705 Double_t rad2 = rxy2+point[2]*point[2];
706 Double_t r=TMath::Sqrt(rad2);
707 Bool_t rzero=kFALSE;
708 if (r<=1E-20) rzero=kTRUE;
709 //localize theta
710 Double_t phi=0;;
711 Double_t th=0.;
712 if (TestShapeBit(kGeoThetaSeg) && (!rzero)) {
713 th = TMath::ACos(point[2]/r)*TMath::RadToDeg();
714 }
715 //localize phi
716 if (TestShapeBit(kGeoPhiSeg)) {
717 phi=TMath::ATan2(point[1], point[0])*TMath::RadToDeg();
718 if (phi<0) phi+=360.;
719 }
720 if (iact<3 && safe) {
721 saf[0]=(TGeoShape::IsSameWithinTolerance(fRmin,0))?TGeoShape::Big():r-fRmin;
722 saf[1]=fRmax-r;
723 saf[2]=saf[3]=saf[4]=saf[5]= TGeoShape::Big();
724 if (TestShapeBit(kGeoThetaSeg)) {
725 if (fTheta1>0) {
726 saf[2] = r*TMath::Sin((th-fTheta1)*TMath::DegToRad());
727 }
728 if (fTheta2<180) {
729 saf[3] = r*TMath::Sin((fTheta2-th)*TMath::DegToRad());
730 }
731 }
732 if (TestShapeBit(kGeoPhiSeg)) {
733 Double_t dph1=phi-fPhi1;
734 if (dph1<0) dph1+=360.;
735 if (dph1<=90.) saf[4]=rxy*TMath::Sin(dph1*TMath::DegToRad());
736 Double_t dph2=fPhi2-phi;
737 if (dph2<0) dph2+=360.;
738 if (dph2<=90.) saf[5]=rxy*TMath::Sin(dph2*TMath::DegToRad());
739 }
740 *safe = saf[TMath::LocMin(6, &saf[0])];
741 if (iact==0) return TGeoShape::Big();
742 if (iact==1 && step<*safe) return TGeoShape::Big();
743 }
744 // compute distance to shape
745 if (rzero) {
746// gGeoManager->SetNormalChecked(1.);
747 return fRmax;
748 }
749 // first do rmin, rmax
750 Double_t b,delta, xnew,ynew,znew, phi0, ddp;
751 Double_t rdotn = point[0]*dir[0]+point[1]*dir[1]+point[2]*dir[2];
752 Double_t sn1 = TGeoShape::Big();
753 // Inner sphere
754 if (fRmin>0) {
755 // Protection in case point is actually outside the sphere
756 if (r <= fRmin+TGeoShape::Tolerance()) {
757 if (rdotn<0) return 0.0;
758 } else {
759 if (rdotn<0) sn1 = DistToSphere(point, dir, fRmin, kFALSE);
760 }
761 }
762 // Outer sphere
763 if (r >= fRmax-TGeoShape::Tolerance()) {
764 if (rdotn>=0) return 0.0;
765 }
766 Double_t sn2 = DistToSphere(point, dir, fRmax, kFALSE, kFALSE);
767 Double_t sr = TMath::Min(sn1, sn2);
768 // check theta conical surfaces
769 sn1 = sn2 = TGeoShape::Big();
770 if (TestShapeBit(kGeoThetaSeg)) {
771 if (TGeoShape::IsSameWithinTolerance(fTheta1,90)) {
772 // surface is a plane
773 if (point[2]*dir[2]<0) sn1 = -point[2]/dir[2];
774 } else {
775 if (fTheta1>0) {
776 Double_t r1,r2,z1,z2,dz,ptnew[3];
777 Double_t si = TMath::Sin(fTheta1*TMath::DegToRad());
778 Double_t ci = TMath::Cos(fTheta1*TMath::DegToRad());
779 if (ci>0) {
780 r1 = fRmin*si;
781 z1 = fRmin*ci;
782 r2 = fRmax*si;
783 z2 = fRmax*ci;
784 } else {
785 r1 = fRmax*si;
786 z1 = fRmax*ci;
787 r2 = fRmin*si;
788 z2 = fRmin*ci;
789 }
790 dz = 0.5*(z2-z1);
791 ptnew[0] = point[0];
792 ptnew[1] = point[1];
793 ptnew[2] = point[2]-0.5*(z1+z2);
794 Double_t zinv = 1./dz;
795 Double_t rin = 0.5*(r1+r2+(r2-r1)*ptnew[2]*zinv);
796 // Protection in case point is outside
797 Double_t sigz = TMath::Sign(1.,point[2]);
798 if (sigz*ci>0 && sigz*rxy2 < sigz*rin*(rin+sigz*TGeoShape::Tolerance())) {
799 Double_t ddotn = ptnew[0]*dir[0]+ptnew[1]*dir[1]+0.5*(r1-r2)*dir[2]*zinv*TMath::Sqrt(rxy2);
800 if (sigz*ddotn<=0) return 0.0;
801 } else {
802 TGeoCone::DistToCone(ptnew, dir, dz, r1, r2, b, delta);
803 if (delta>0) {
804 Double_t snxt = -b-delta;
805 znew = ptnew[2]+snxt*dir[2];
806 if (snxt>0 && TMath::Abs(znew)<dz) {
807 if (!TestShapeBit(kGeoPhiSeg)) sn1=snxt;
808 else {
809 xnew = ptnew[0]+snxt*dir[0];
810 ynew = ptnew[1]+snxt*dir[1];
811 phi0 = TMath::ATan2(ynew, xnew) * TMath::RadToDeg();
812 ddp = phi0-fPhi1;
813 while (ddp<0) ddp+=360.;
814 if (ddp<=fPhi2-fPhi1) sn1 = snxt;
815 }
816 }
817 if (sn1>1E10) {
818 snxt = -b+delta;
819 znew = ptnew[2]+snxt*dir[2];
820 if (snxt>0 && TMath::Abs(znew)<dz) {
821 if (!TestShapeBit(kGeoPhiSeg)) sn1=snxt;
822 else {
823 xnew = ptnew[0]+snxt*dir[0];
824 ynew = ptnew[1]+snxt*dir[1];
825 phi0 = TMath::ATan2(ynew, xnew) * TMath::RadToDeg();
826 ddp = phi0-fPhi1;
827 while (ddp<0) ddp+=360.;
828 if (ddp<=fPhi2-fPhi1) sn1 = snxt;
829 }
830 }
831 }
832 }
833 }
834 }
835 }
836 if (TGeoShape::IsSameWithinTolerance(fTheta2,90)) {
837 // surface is a plane
838 if (point[2]*dir[2]<0) sn1 = -point[2]/dir[2];
839 } else {
840 if (fTheta2<180) {
841 Double_t r1,r2,z1,z2,dz,ptnew[3];
842 Double_t si = TMath::Sin(fTheta2*TMath::DegToRad());
843 Double_t ci = TMath::Cos(fTheta2*TMath::DegToRad());
844 if (ci>0) {
845 r1 = fRmin*si;
846 z1 = fRmin*ci;
847 r2 = fRmax*si;
848 z2 = fRmax*ci;
849 } else {
850 r1 = fRmax*si;
851 z1 = fRmax*ci;
852 r2 = fRmin*si;
853 z2 = fRmin*ci;
854 }
855 dz = 0.5*(z2-z1);
856 ptnew[0] = point[0];
857 ptnew[1] = point[1];
858 ptnew[2] = point[2]-0.5*(z1+z2);
859 Double_t zinv = 1./dz;
860 Double_t rin = 0.5*(r1+r2+(r2-r1)*ptnew[2]*zinv);
861 // Protection in case point is outside
862 Double_t sigz = TMath::Sign(1.,point[2]);
863 if (sigz*ci>0 && sigz*rxy2 > sigz*rin*(rin-sigz*TGeoShape::Tolerance())) {
864 Double_t ddotn = ptnew[0]*dir[0]+ptnew[1]*dir[1]+0.5*(r1-r2)*dir[2]*zinv*TMath::Sqrt(rxy2);
865 if (sigz*ddotn>=0) return 0.0;
866 } else {
867 TGeoCone::DistToCone(ptnew, dir, dz, r1, r2, b, delta);
868 if (delta>0) {
869 Double_t snxt = -b-delta;
870 znew = ptnew[2]+snxt*dir[2];
871 if (snxt>0 && TMath::Abs(znew)<dz) {
872 if (!TestShapeBit(kGeoPhiSeg)) sn2 = snxt;
873 else {
874 xnew = ptnew[0]+snxt*dir[0];
875 ynew = ptnew[1]+snxt*dir[1];
876 phi0 = TMath::ATan2(ynew, xnew) * TMath::RadToDeg();
877 ddp = phi0-fPhi1;
878 while (ddp<0) ddp+=360.;
879 if (ddp<=fPhi2-fPhi1) sn2 = snxt;
880 }
881 }
882 if (sn2>1E10) {
883 snxt = -b+delta;
884 znew = ptnew[2]+snxt*dir[2];
885 if (snxt>0 && TMath::Abs(znew)<dz) {
886 if (!TestShapeBit(kGeoPhiSeg)) sn2 = snxt;
887 else {
888 xnew = ptnew[0]+snxt*dir[0];
889 ynew = ptnew[1]+snxt*dir[1];
890 phi0 = TMath::ATan2(ynew, xnew) * TMath::RadToDeg();
891 ddp = phi0-fPhi1;
892 while (ddp<0) ddp+=360.;
893 if (ddp<=fPhi2-fPhi1) sn2 = snxt;
894 }
895 }
896 }
897 }
898 }
899 }
900 }
901 }
902 Double_t st = TMath::Min(sn1,sn2);
903 Double_t sp = TGeoShape::Big();
904 if (TestShapeBit(kGeoPhiSeg)) {
905 Double_t s1 = TMath::Sin(fPhi1*TMath::DegToRad());
906 Double_t c1 = TMath::Cos(fPhi1*TMath::DegToRad());
907 Double_t s2 = TMath::Sin(fPhi2*TMath::DegToRad());
908 Double_t c2 = TMath::Cos(fPhi2*TMath::DegToRad());
909 Double_t phim = 0.5*(fPhi1+fPhi2);
910 Double_t sm = TMath::Sin(phim*TMath::DegToRad());
911 Double_t cm = TMath::Cos(phim*TMath::DegToRad());
912 sp = TGeoShape::DistToPhiMin(point, dir, s1, c1, s2, c2, sm, cm);
913 }
914 Double_t snxt = TMath::Min(sr, st);
915 snxt = TMath::Min(snxt, sp);
916 return snxt;
917}
918
919////////////////////////////////////////////////////////////////////////////////
920/// compute distance to sphere of radius rsph. Direction has to be a unit vector
921
922Double_t TGeoSphere::DistToSphere(const Double_t *point, const Double_t *dir, Double_t rsph, Bool_t check, Bool_t firstcross) const
923{
924 if (rsph<=0) return 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 Double_t s;
937 if (in) {
938 s=-b+d;
939 } else {
940 s = (firstcross)?(-b-d):(-b+d);
941 }
942 if (s<0) return TGeoShape::Big();
943 if (!check) return s;
944 for (i=0; i<3; i++) pt[i]=point[i]+s*dir[i];
945 // check theta and phi ranges
946 if (IsPointInside(&pt[0], kFALSE)) return s;
947 return TGeoShape::Big();
948}
949
950////////////////////////////////////////////////////////////////////////////////
951
952TGeoVolume *TGeoSphere::Divide(TGeoVolume * voldiv, const char * divname, Int_t iaxis, Int_t ndiv,
953 Double_t start, Double_t step)
954{
955 TGeoShape *shape; //--- shape to be created
956 TGeoVolume *vol; //--- division volume to be created
957 TGeoVolumeMulti *vmulti; //--- generic divided volume
958 TGeoPatternFinder *finder; //--- finder to be attached
959 TString opt = ""; //--- option to be attached
960 Int_t id;
961 Double_t end = start+ndiv*step;
962 switch (iaxis) {
963 case 1: //--- R division
964 finder = new TGeoPatternSphR(voldiv, ndiv, start, end);
965 vmulti = gGeoManager->MakeVolumeMulti(divname, voldiv->GetMedium());
966 voldiv->SetFinder(finder);
967 finder->SetDivIndex(voldiv->GetNdaughters());
968 for (id=0; id<ndiv; id++) {
969 shape = new TGeoSphere(start+id*step, start+(id+1)*step, fTheta1,fTheta2,fPhi1,fPhi2);
970 vol = new TGeoVolume(divname, shape, voldiv->GetMedium());
971 vmulti->AddVolume(vol);
972 opt = "R";
973 voldiv->AddNodeOffset(vol, id, 0, opt.Data());
974 ((TGeoNodeOffset*)voldiv->GetNodes()->At(voldiv->GetNdaughters()-1))->SetFinder(finder);
975 }
976 return vmulti;
977 case 2: //--- Phi division
978 finder = new TGeoPatternSphPhi(voldiv, ndiv, start, end);
979 voldiv->SetFinder(finder);
980 finder->SetDivIndex(voldiv->GetNdaughters());
981 shape = new TGeoSphere(fRmin, fRmax, fTheta1, fTheta2, -step/2, step/2);
982 vol = new TGeoVolume(divname, shape, voldiv->GetMedium());
983 vmulti = gGeoManager->MakeVolumeMulti(divname, voldiv->GetMedium());
984 vmulti->AddVolume(vol);
985 opt = "Phi";
986 for (id=0; id<ndiv; id++) {
987 voldiv->AddNodeOffset(vol, id, start+id*step+step/2, opt.Data());
988 ((TGeoNodeOffset*)voldiv->GetNodes()->At(voldiv->GetNdaughters()-1))->SetFinder(finder);
989 }
990 return vmulti;
991 case 3: //--- Theta division
992 finder = new TGeoPatternSphTheta(voldiv, ndiv, start, end);
993 vmulti = gGeoManager->MakeVolumeMulti(divname, voldiv->GetMedium());
994 voldiv->SetFinder(finder);
995 finder->SetDivIndex(voldiv->GetNdaughters());
996 for (id=0; id<ndiv; id++) {
997 shape = new TGeoSphere(fRmin,fRmax,start+id*step, start+(id+1)*step,fPhi1,fPhi2);
998 vol = new TGeoVolume(divname, shape, voldiv->GetMedium());
999 vmulti->AddVolume(vol);
1000 opt = "Theta";
1001 voldiv->AddNodeOffset(vol, id, 0, opt.Data());
1002 ((TGeoNodeOffset*)voldiv->GetNodes()->At(voldiv->GetNdaughters()-1))->SetFinder(finder);
1003 }
1004 return vmulti;
1005 default:
1006 Error("Divide", "In shape %s wrong axis type for division", GetName());
1007 return 0;
1008 }
1009}
1010
1011////////////////////////////////////////////////////////////////////////////////
1012/// Returns name of axis IAXIS.
1013
1014const char *TGeoSphere::GetAxisName(Int_t iaxis) const
1015{
1016 switch (iaxis) {
1017 case 1:
1018 return "R";
1019 case 2:
1020 return "PHI";
1021 case 3:
1022 return "THETA";
1023 default:
1024 return "UNDEFINED";
1025 }
1026}
1027
1028////////////////////////////////////////////////////////////////////////////////
1029/// Get range of shape for a given axis.
1030
1031Double_t TGeoSphere::GetAxisRange(Int_t iaxis, Double_t &xlo, Double_t &xhi) const
1032{
1033 xlo = 0;
1034 xhi = 0;
1035 Double_t dx = 0;
1036 switch (iaxis) {
1037 case 1:
1038 xlo = fRmin;
1039 xhi = fRmax;
1040 dx = xhi-xlo;
1041 return dx;
1042 case 2:
1043 xlo = fPhi1;
1044 xhi = fPhi2;
1045 dx = xhi-xlo;
1046 return dx;
1047 case 3:
1048 xlo = fTheta1;
1049 xhi = fTheta2;
1050 dx = xhi-xlo;
1051 return dx;
1052 }
1053 return dx;
1054}
1055
1056////////////////////////////////////////////////////////////////////////////////
1057/// Fill vector param[4] with the bounding cylinder parameters. The order
1058/// is the following : Rmin, Rmax, Phi1, Phi2
1059
1060void TGeoSphere::GetBoundingCylinder(Double_t *param) const
1061{
1062 Double_t smin = TMath::Sin(fTheta1*TMath::DegToRad());
1063 Double_t smax = TMath::Sin(fTheta2*TMath::DegToRad());
1064 if (smin>smax) {
1065 Double_t a = smin;
1066 smin = smax;
1067 smax = a;
1068 }
1069 param[0] = fRmin*smin; // Rmin
1070 param[0] *= param[0];
1071 if (((90.-fTheta1)*(fTheta2-90.))>=0) smax = 1.;
1072 param[1] = fRmax*smax; // Rmax
1073 param[1] *= param[1];
1074 param[2] = (fPhi1<0)?(fPhi1+360.):fPhi1; // Phi1
1075 param[3] = fPhi2;
1076 if (TGeoShape::IsSameWithinTolerance(param[3]-param[2],360)) { // Phi2
1077 param[2] = 0.;
1078 param[3] = 360.;
1079 }
1080 while (param[3]<param[2]) param[3]+=360.;
1081}
1082
1083////////////////////////////////////////////////////////////////////////////////
1084/// print shape parameters
1085
1086void TGeoSphere::InspectShape() const
1087{
1088 printf("*** Shape %s: TGeoSphere ***\n", GetName());
1089 printf(" Rmin = %11.5f\n", fRmin);
1090 printf(" Rmax = %11.5f\n", fRmax);
1091 printf(" Th1 = %11.5f\n", fTheta1);
1092 printf(" Th2 = %11.5f\n", fTheta2);
1093 printf(" Ph1 = %11.5f\n", fPhi1);
1094 printf(" Ph2 = %11.5f\n", fPhi2);
1095 printf(" Bounding box:\n");
1096 TGeoBBox::InspectShape();
1097}
1098
1099////////////////////////////////////////////////////////////////////////////////
1100/// Creates a TBuffer3D describing *this* shape.
1101/// Coordinates are in local reference frame.
1102
1103TBuffer3D *TGeoSphere::MakeBuffer3D() const
1104{
1105 Bool_t full = kTRUE;
1106 if (TestShapeBit(kGeoThetaSeg) || TestShapeBit(kGeoPhiSeg)) full = kFALSE;
1107 Int_t ncenter = 1;
1108 if (full || TestShapeBit(kGeoRSeg)) ncenter = 0;
1109 Int_t nup = (fTheta1>0)?0:1;
1110 Int_t ndown = (fTheta2<180)?0:1;
1111 // number of different latitudes, excluding 0 and 180 degrees
1112 Int_t nlat = fNz+1-(nup+ndown);
1113 // number of different longitudes
1114 Int_t nlong = fNseg;
1115 if (TestShapeBit(kGeoPhiSeg)) nlong++;
1116
1117 Int_t nbPnts = nlat*nlong+nup+ndown+ncenter;
1118 if (TestShapeBit(kGeoRSeg)) nbPnts *= 2;
1119
1120 Int_t nbSegs = nlat*fNseg + (nlat-1+nup+ndown)*nlong; // outer sphere
1121 if (TestShapeBit(kGeoRSeg)) nbSegs *= 2; // inner sphere
1122 if (TestShapeBit(kGeoPhiSeg)) nbSegs += 2*nlat+nup+ndown; // 2 phi planes
1123 nbSegs += nlong * (2-nup - ndown); // connecting cones
1124
1125 Int_t nbPols = fNz*fNseg; // outer
1126 if (TestShapeBit(kGeoRSeg)) nbPols *=2; // inner
1127 if (TestShapeBit(kGeoPhiSeg)) nbPols += 2*fNz; // 2 phi planes
1128 nbPols += (2-nup-ndown)*fNseg; // connecting
1129
1130 TBuffer3D* buff = new TBuffer3D(TBuffer3DTypes::kGeneric,
1131 nbPnts, 3*nbPnts, nbSegs, 3*nbSegs, nbPols, 6*nbPols);
1132
1133 if (buff)
1134 {
1135 SetPoints(buff->fPnts);
1136 SetSegsAndPols(*buff);
1137 }
1138
1139 return buff;
1140}
1141
1142////////////////////////////////////////////////////////////////////////////////
1143/// Fill TBuffer3D structure for segments and polygons.
1144
1145void TGeoSphere::SetSegsAndPols(TBuffer3D & buff) const
1146{
1147 // Bool_t full = kTRUE;
1148 // if (TestShapeBit(kGeoThetaSeg) || TestShapeBit(kGeoPhiSeg)) full = kFALSE;
1149 // Int_t ncenter = 1;
1150 // if (full || TestShapeBit(kGeoRSeg)) ncenter = 0;
1151 Int_t nup = (fTheta1 > 0) ? 0 : 1;
1152 Int_t ndown = (fTheta2 < 180) ? 0 : 1;
1153 // number of different latitudes, excluding 0 and 180 degrees
1154 Int_t nlat = fNz+1-(nup+ndown);
1155 // number of different longitudes
1156 Int_t nlong = fNseg;
1157 if (TestShapeBit(kGeoPhiSeg)) nlong++;
1158
1159 // Int_t nbPnts = nlat*nlong+nup+ndown+ncenter;
1160 // if (TestShapeBit(kGeoRSeg)) nbPnts *= 2;
1161
1162 // Int_t nbSegs = nlat*fNseg + (nlat-1+nup+ndown)*nlong; // outer sphere
1163 // if (TestShapeBit(kGeoRSeg)) nbSegs *= 2; // inner sphere
1164 // if (TestShapeBit(kGeoPhiSeg)) nbSegs += 2*nlat+nup+ndown; // 2 phi planes
1165 // nbSegs += nlong * (2-nup - ndown); // connecting cones
1166
1167 // Int_t nbPols = fNz*fNseg; // outer
1168 // if (TestShapeBit(kGeoRSeg)) nbPols *=2; // inner
1169 // if (TestShapeBit(kGeoPhiSeg)) nbPols += 2*fNz; // 2 phi planes
1170 // nbPols += (2-nup-ndown)*fNseg; // connecting
1171
1172 Int_t c = GetBasicColor();
1173 Int_t i, j;
1174 Int_t indx;
1175 indx = 0;
1176 // outside sphere
1177 // loop all segments on latitudes (except 0 and 180 degrees)
1178 // [0, nlat*fNseg)
1179 Int_t indpar = 0;
1180 for (i=0; i<nlat; i++) {
1181 for (j=0; j<fNseg; j++) {
1182 buff.fSegs[indx++] = c;
1183 buff.fSegs[indx++] = i*nlong+j;
1184 buff.fSegs[indx++] = i*nlong+(j+1)%nlong;
1185 }
1186 }
1187 // loop all segments on longitudes
1188 // nlat*fNseg + [0, (nlat-1)*nlong)
1189 Int_t indlong = indpar + nlat*fNseg;
1190 for (i=0; i<nlat-1; i++) {
1191 for (j=0; j<nlong; j++) {
1192 buff.fSegs[indx++] = c;
1193 buff.fSegs[indx++] = i*nlong+j;
1194 buff.fSegs[indx++] = (i+1)*nlong+j;
1195 }
1196 }
1197 Int_t indup = indlong + (nlat-1)*nlong;
1198 // extra longitudes on top
1199 // nlat*fNseg+(nlat-1)*nlong + [0, nlong)
1200 if (nup) {
1201 Int_t indpup = nlat*nlong;
1202 for (j=0; j<nlong; j++) {
1203 buff.fSegs[indx++] = c;
1204 buff.fSegs[indx++] = j;
1205 buff.fSegs[indx++] = indpup;
1206 }
1207 }
1208 Int_t inddown = indup + nup*nlong;
1209 // extra longitudes on bottom
1210 // nlat*fNseg+(nlat+nup-1)*nlong + [0, nlong)
1211 if (ndown) {
1212 Int_t indpdown = nlat*nlong+nup;
1213 for (j=0; j<nlong; j++) {
1214 buff.fSegs[indx++] = c;
1215 buff.fSegs[indx++] = (nlat-1)*nlong+j;
1216 buff.fSegs[indx++] = indpdown;
1217 }
1218 }
1219 Int_t indparin = inddown + ndown*nlong;
1220 Int_t indlongin = indparin;
1221 Int_t indupin = indparin;
1222 Int_t inddownin = indparin;
1223 Int_t indphi = indparin;
1224 // inner sphere
1225 Int_t indptin = nlat*nlong + nup + ndown;
1226 Int_t iptcenter = indptin;
1227 // nlat*fNseg+(nlat+nup+ndown-1)*nlong
1228 if (TestShapeBit(kGeoRSeg)) {
1229 indlongin = indparin + nlat*fNseg;
1230 indupin = indlongin + (nlat-1)*nlong;
1231 inddownin = indupin + nup*nlong;
1232 // loop all segments on latitudes (except 0 and 180 degrees)
1233 // indsegin + [0, nlat*fNseg)
1234 for (i=0; i<nlat; i++) {
1235 for (j=0; j<fNseg; j++) {
1236 buff.fSegs[indx++] = c+1;
1237 buff.fSegs[indx++] = indptin + i*nlong+j;
1238 buff.fSegs[indx++] = indptin + i*nlong+(j+1)%nlong;
1239 }
1240 }
1241 // loop all segments on longitudes
1242 // indsegin + nlat*fNseg + [0, (nlat-1)*nlong)
1243 for (i=0; i<nlat-1; i++) {
1244 for (j=0; j<nlong; j++) {
1245 buff.fSegs[indx++] = c+1;
1246 buff.fSegs[indx++] = indptin + i*nlong+j;
1247 buff.fSegs[indx++] = indptin + (i+1)*nlong+j;
1248 }
1249 }
1250 // extra longitudes on top
1251 // indsegin + nlat*fNseg+(nlat-1)*nlong + [0, nlong)
1252 if (nup) {
1253 Int_t indupltop = indptin + nlat*nlong;
1254 for (j=0; j<nlong; j++) {
1255 buff.fSegs[indx++] = c+1;
1256 buff.fSegs[indx++] = indptin + j;
1257 buff.fSegs[indx++] = indupltop;
1258 }
1259 }
1260 // extra longitudes on bottom
1261 // indsegin + nlat*fNseg+(nlat+nup-1)*nlong + [0, nlong)
1262 if (ndown) {
1263 Int_t indpdown = indptin + nlat*nlong+nup;
1264 for (j=0; j<nlong; j++) {
1265 buff.fSegs[indx++] = c+1;
1266 buff.fSegs[indx++] = indptin + (nlat-1)*nlong+j;
1267 buff.fSegs[indx++] = indpdown;
1268 }
1269 }
1270 indphi = inddownin + ndown*nlong;
1271 }
1272 Int_t indtheta = indphi;
1273 // Segments on phi planes
1274 if (TestShapeBit(kGeoPhiSeg)) {
1275 indtheta += 2*nlat + nup + ndown;
1276 for (j=0; j<nlat; j++) {
1277 buff.fSegs[indx++] = c+2;
1278 buff.fSegs[indx++] = j*nlong;
1279 if (TestShapeBit(kGeoRSeg)) buff.fSegs[indx++] = indptin + j*nlong;
1280 else buff.fSegs[indx++] = iptcenter;
1281 }
1282 for (j=0; j<nlat; j++) {
1283 buff.fSegs[indx++] = c+2;
1284 buff.fSegs[indx++] = (j+1)*nlong-1;
1285 if (TestShapeBit(kGeoRSeg)) buff.fSegs[indx++] = indptin + (j+1)*nlong-1;
1286 else buff.fSegs[indx++] = iptcenter;
1287 }
1288 if (nup) {
1289 buff.fSegs[indx++] = c+2;
1290 buff.fSegs[indx++] = nlat*nlong;
1291 if (TestShapeBit(kGeoRSeg)) buff.fSegs[indx++] = indptin + nlat*nlong;
1292 else buff.fSegs[indx++] = iptcenter;
1293 }
1294 if (ndown) {
1295 buff.fSegs[indx++] = c+2;
1296 buff.fSegs[indx++] = nlat*nlong+nup;
1297 if (TestShapeBit(kGeoRSeg)) buff.fSegs[indx++] = indptin + nlat*nlong+nup;
1298 else buff.fSegs[indx++] = iptcenter;
1299 }
1300 }
1301 // Segments on cones
1302 if (!nup) {
1303 for (j=0; j<nlong; j++) {
1304 buff.fSegs[indx++] = c+2;
1305 buff.fSegs[indx++] = j;
1306 if (TestShapeBit(kGeoRSeg)) buff.fSegs[indx++] = indptin + j;
1307 else buff.fSegs[indx++] = iptcenter;
1308 }
1309 }
1310 if (!ndown) {
1311 for (j=0; j<nlong; j++) {
1312 buff.fSegs[indx++] = c+2;
1313 buff.fSegs[indx++] = (nlat-1)*nlong + j;
1314 if (TestShapeBit(kGeoRSeg)) buff.fSegs[indx++] = indptin + (nlat-1)*nlong +j;
1315 else buff.fSegs[indx++] = iptcenter;
1316 }
1317 }
1318
1319 indx = 0;
1320 // Fill polygons for outside sphere (except 0/180)
1321 for (i=0; i<nlat-1; i++) {
1322 for (j=0; j<fNseg; j++) {
1323 buff.fPols[indx++] = c;
1324 buff.fPols[indx++] = 4;
1325 buff.fPols[indx++] = indpar+i*fNseg+j;
1326 buff.fPols[indx++] = indlong+i*nlong+(j+1)%nlong;
1327 buff.fPols[indx++] = indpar+(i+1)*fNseg+j;
1328 buff.fPols[indx++] = indlong+i*nlong+j;
1329 }
1330 }
1331 // upper
1332 if (nup) {
1333 for (j=0; j<fNseg; j++) {
1334 buff.fPols[indx++] = c;
1335 buff.fPols[indx++] = 3;
1336 buff.fPols[indx++] = indup + j;
1337 buff.fPols[indx++] = indup + (j+1)%nlong;
1338 buff.fPols[indx++] = indpar + j;
1339 }
1340 }
1341 // lower
1342 if (ndown) {
1343 for (j=0; j<fNseg; j++) {
1344 buff.fPols[indx++] = c;
1345 buff.fPols[indx++] = 3;
1346 buff.fPols[indx++] = inddown + j;
1347 buff.fPols[indx++] = indpar + (nlat-1)*fNseg + j;
1348 buff.fPols[indx++] = inddown + (j+1)%nlong;
1349 }
1350 }
1351 // Fill polygons for inside sphere (except 0/180)
1352
1353 if (TestShapeBit(kGeoRSeg)) {
1354 for (i=0; i<nlat-1; i++) {
1355 for (j=0; j<fNseg; j++) {
1356 buff.fPols[indx++] = c+1;
1357 buff.fPols[indx++] = 4;
1358 buff.fPols[indx++] = indparin+i*fNseg+j;
1359 buff.fPols[indx++] = indlongin+i*nlong+j;
1360 buff.fPols[indx++] = indparin+(i+1)*fNseg+j;
1361 buff.fPols[indx++] = indlongin+i*nlong+(j+1)%nlong;
1362 }
1363 }
1364 // upper
1365 if (nup) {
1366 for (j=0; j<fNseg; j++) {
1367 buff.fPols[indx++] = c+1;
1368 buff.fPols[indx++] = 3;
1369 buff.fPols[indx++] = indupin + j;
1370 buff.fPols[indx++] = indparin + j;
1371 buff.fPols[indx++] = indupin + (j+1)%nlong;
1372 }
1373 }
1374 // lower
1375 if (ndown) {
1376 for (j=0; j<fNseg; j++) {
1377 buff.fPols[indx++] = c+1;
1378 buff.fPols[indx++] = 3;
1379 buff.fPols[indx++] = inddownin + j;
1380 buff.fPols[indx++] = inddownin + (j+1)%nlong;
1381 buff.fPols[indx++] = indparin + (nlat-1)*fNseg + j;
1382 }
1383 }
1384 }
1385 // Polygons on phi planes
1386 if (TestShapeBit(kGeoPhiSeg)) {
1387 for (i=0; i<nlat-1; i++) {
1388 buff.fPols[indx++] = c+2;
1389 if (TestShapeBit(kGeoRSeg)) {
1390 buff.fPols[indx++] = 4;
1391 buff.fPols[indx++] = indlong + i*nlong;
1392 buff.fPols[indx++] = indphi + i + 1;
1393 buff.fPols[indx++] = indlongin + i*nlong;
1394 buff.fPols[indx++] = indphi + i;
1395 } else {
1396 buff.fPols[indx++] = 3;
1397 buff.fPols[indx++] = indlong + i*nlong;
1398 buff.fPols[indx++] = indphi + i + 1;
1399 buff.fPols[indx++] = indphi + i;
1400 }
1401 }
1402 for (i=0; i<nlat-1; i++) {
1403 buff.fPols[indx++] = c+2;
1404 if (TestShapeBit(kGeoRSeg)) {
1405 buff.fPols[indx++] = 4;
1406 buff.fPols[indx++] = indlong + (i+1)*nlong-1;
1407 buff.fPols[indx++] = indphi + nlat + i;
1408 buff.fPols[indx++] = indlongin + (i+1)*nlong-1;
1409 buff.fPols[indx++] = indphi + nlat + i + 1;
1410 } else {
1411 buff.fPols[indx++] = 3;
1412 buff.fPols[indx++] = indlong + (i+1)*nlong-1;
1413 buff.fPols[indx++] = indphi + nlat + i;
1414 buff.fPols[indx++] = indphi + nlat + i + 1;
1415 }
1416 }
1417 if (nup) {
1418 buff.fPols[indx++] = c+2;
1419 if (TestShapeBit(kGeoRSeg)) {
1420 buff.fPols[indx++] = 4;
1421 buff.fPols[indx++] = indup;
1422 buff.fPols[indx++] = indphi;
1423 buff.fPols[indx++] = indupin;
1424 buff.fPols[indx++] = indphi + 2*nlat;
1425 } else {
1426 buff.fPols[indx++] = 3;
1427 buff.fPols[indx++] = indup;
1428 buff.fPols[indx++] = indphi;
1429 buff.fPols[indx++] = indphi + 2*nlat;
1430 }
1431 buff.fPols[indx++] = c+2;
1432 if (TestShapeBit(kGeoRSeg)) {
1433 buff.fPols[indx++] = 4;
1434 buff.fPols[indx++] = indup+nlong-1;
1435 buff.fPols[indx++] = indphi + 2*nlat;
1436 buff.fPols[indx++] = indupin+nlong-1;
1437 buff.fPols[indx++] = indphi + nlat;
1438 } else {
1439 buff.fPols[indx++] = 3;
1440 buff.fPols[indx++] = indup+nlong-1;
1441 buff.fPols[indx++] = indphi + 2*nlat;
1442 buff.fPols[indx++] = indphi + nlat;
1443 }
1444 }
1445 if (ndown) {
1446 buff.fPols[indx++] = c+2;
1447 if (TestShapeBit(kGeoRSeg)) {
1448 buff.fPols[indx++] = 4;
1449 buff.fPols[indx++] = inddown;
1450 buff.fPols[indx++] = indphi + 2*nlat + nup;
1451 buff.fPols[indx++] = inddownin;
1452 buff.fPols[indx++] = indphi + nlat-1;
1453 } else {
1454 buff.fPols[indx++] = 3;
1455 buff.fPols[indx++] = inddown;
1456 buff.fPols[indx++] = indphi + 2*nlat + nup;
1457 buff.fPols[indx++] = indphi + nlat-1;
1458 }
1459 buff.fPols[indx++] = c+2;
1460 if (TestShapeBit(kGeoRSeg)) {
1461 buff.fPols[indx++] = 4;
1462 buff.fPols[indx++] = inddown+nlong-1;
1463 buff.fPols[indx++] = indphi + 2*nlat-1;
1464 buff.fPols[indx++] = inddownin+nlong-1;
1465 buff.fPols[indx++] = indphi + 2*nlat+nup;
1466 } else {
1467 buff.fPols[indx++] = 3;
1468 buff.fPols[indx++] = inddown+nlong-1;
1469 buff.fPols[indx++] = indphi + 2*nlat-1;
1470 buff.fPols[indx++] = indphi + 2*nlat+nup;
1471 }
1472 }
1473 }
1474 // Polygons on cones
1475 if (!nup) {
1476 for (j=0; j<fNseg; j++) {
1477 buff.fPols[indx++] = c+2;
1478 if (TestShapeBit(kGeoRSeg)) {
1479 buff.fPols[indx++] = 4;
1480 buff.fPols[indx++] = indpar+j;
1481 buff.fPols[indx++] = indtheta + j;
1482 buff.fPols[indx++] = indparin + j;
1483 buff.fPols[indx++] = indtheta + (j+1)%nlong;
1484 } else {
1485 buff.fPols[indx++] = 3;
1486 buff.fPols[indx++] = indpar+j;
1487 buff.fPols[indx++] = indtheta + j;
1488 buff.fPols[indx++] = indtheta + (j+1)%nlong;
1489 }
1490 }
1491 }
1492 if (!ndown) {
1493 for (j=0; j<fNseg; j++) {
1494 buff.fPols[indx++] = c+2;
1495 if (TestShapeBit(kGeoRSeg)) {
1496 buff.fPols[indx++] = 4;
1497 buff.fPols[indx++] = indpar+(nlat-1)*fNseg+j;
1498 buff.fPols[indx++] = indtheta + (1-nup)*nlong +(j+1)%nlong;
1499 buff.fPols[indx++] = indparin + (nlat-1)*fNseg + j;
1500 buff.fPols[indx++] = indtheta + (1-nup)*nlong + j;
1501 } else {
1502 buff.fPols[indx++] = 3;
1503 buff.fPols[indx++] = indpar+(nlat-1)*fNseg+j;
1504 buff.fPols[indx++] = indtheta + (1-nup)*nlong +(j+1)%nlong;
1505 buff.fPols[indx++] = indtheta + (1-nup)*nlong + j;
1506 }
1507 }
1508 }
1509}
1510
1511////////////////////////////////////////////////////////////////////////////////
1512/// computes the closest distance from given point to this shape, according
1513/// to option. The matching point on the shape is stored in spoint.
1514
1515Double_t TGeoSphere::Safety(const Double_t *point, Bool_t in) const
1516{
1517 Double_t r2 = point[0]*point[0]+point[1]*point[1]+point[2]*point[2];
1518 Double_t r=TMath::Sqrt(r2);
1519 Bool_t rzero=kFALSE;
1520 if (r<=1E-20) rzero=kTRUE;
1521 //localize theta
1522 Double_t th=0.;
1523 if (TestShapeBit(kGeoThetaSeg) && (!rzero)) {
1524 th = TMath::ACos(point[2]/r)*TMath::RadToDeg();
1525 }
1526 Double_t saf[4];
1527 saf[0]=(TGeoShape::IsSameWithinTolerance(fRmin,0) && !TestShapeBit(kGeoThetaSeg) && !TestShapeBit(kGeoPhiSeg))?TGeoShape::Big():r-fRmin;
1528 saf[1]=fRmax-r;
1529 saf[2]=saf[3]= TGeoShape::Big();
1530 if (TestShapeBit(kGeoThetaSeg)) {
1531 if (fTheta1>0) saf[2] = r*TMath::Sin((th-fTheta1)*TMath::DegToRad());
1532 if (fTheta2<180) saf[3] = r*TMath::Sin((fTheta2-th)*TMath::DegToRad());
1533 }
1534 Double_t safphi = TGeoShape::Big();
1535 if (TestShapeBit(kGeoPhiSeg)) safphi = TGeoShape::SafetyPhi(point,in,fPhi1,fPhi2);
1536 if (in) {
1537 Double_t 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 Double_t 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
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:101
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:100
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:458
TGeoBBox::InspectShape
virtual void InspectShape() const
Prints shape parameters.
Definition TGeoBBox.cxx:819
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:1059
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:926
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:527
TGeoManager::MakeVolumeMulti
TGeoVolumeMulti * MakeVolumeMulti(const char *name, TGeoMedium *medium)
Make a TGeoVolumeMulti handling a list of volumes.
Definition TGeoManager.cxx:3320
TGeoManager::GetNsegments
Int_t GetNsegments() const
Get number of segments approximating circles.
Definition TGeoManager.cxx:3361
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
TGeoSphere are not just balls having internal and external radii, but sectors of a sphere having defi...
Definition TGeoSphere.h:18
TGeoSphere::~TGeoSphere
virtual ~TGeoSphere()
destructor
Definition TGeoSphere.cxx:127
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:922
TGeoSphere::MakeBuffer3D
virtual TBuffer3D * MakeBuffer3D() const
Creates a TBuffer3D describing this shape.
Definition TGeoSphere.cxx:1103
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:134
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:400
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:1086
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:224
TGeoSphere::ComputeBBox
virtual void ComputeBBox()
compute bounding box of the sphere
Definition TGeoSphere.cxx:149
TGeoSphere::fRmin
Double_t fRmin
Definition TGeoSphere.h:22
TGeoSphere::TGeoSphere
TGeoSphere()
Default constructor.
Definition TGeoSphere.cxx:67
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:332
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:1031
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:1145
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:700
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:952
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:1515
TGeoSphere::GetAxisName
virtual const char * GetAxisName(Int_t iaxis) const
Returns name of axis IAXIS.
Definition TGeoSphere.cxx:1014
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:360
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:287
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:388
TGeoSphere::GetBoundingCylinder
virtual void GetBoundingCylinder(Double_t *param) const
Fill vector param[4] with the bounding cylinder parameters.
Definition TGeoSphere.cxx:1060
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:164
TObject::TestBit
R__ALWAYS_INLINE Bool_t TestBit(UInt_t f) const
Definition TObject.h:201
TObject::ClassName
virtual const char * ClassName() const
Returns name of class to which the object belongs.
Definition TObject.cxx:200
TObject::SetBit
void SetBit(UInt_t f, Bool_t set)
Set or unset the user status bits as specified in f.
Definition TObject.cxx:766
TObject::Error
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition TObject.cxx:963
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:619
TMath::LocMin
Long64_t LocMin(Long64_t n, const T *a)
Return index of array with the minimum element.
Definition TMath.h:922
TMath::Max
Short_t Max(Short_t a, Short_t b)
Definition TMathBase.h:208
TMath::Sign
T1 Sign(T1 a, T2 b)
Definition TMathBase.h:161
TMath::ATan2
Double_t ATan2(Double_t y, Double_t x)
Definition TMath.h:629
TMath::LocMax
Long64_t LocMax(Long64_t n, const T *a)
Return index of array with the maximum element.
Definition TMath.h:950
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:641
TMath::Min
Short_t Min(Short_t a, Short_t b)
Definition TMathBase.h:176
TMath::Cos
Double_t Cos(Double_t)
Definition TMath.h:593
TMath::Pi
constexpr Double_t Pi()
Definition TMath.h:37
TMath::Sin
Double_t Sin(Double_t)
Definition TMath.h:589
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:1168