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