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