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