Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
TGeoTorus.cxx
Go to the documentation of this file.
1 // @(#)root/geom:$Id$
2 // Author: Andrei Gheata 28/07/03
3 
4 /*************************************************************************
5  * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers. *
6  * All rights reserved. *
7  * *
8  * For the licensing terms see $ROOTSYS/LICENSE. *
9  * For the list of contributors see $ROOTSYS/README/CREDITS. *
10  *************************************************************************/
11 
12 /** \class TGeoTorus
13 \ingroup Geometry_classes
14 
15 Torus segment class. A torus has 5 parameters :
16  - R - axial radius
17  - Rmin - inner radius
18  - Rmax - outer radius
19  - Phi1 - starting phi
20  - Dphi - phi extent
21 
22 
23 Begin_Macro(source)
24 {
25  TCanvas *c = new TCanvas("c", "c",0,0,600,600);
26  new TGeoManager("torus", "poza2");
27  TGeoMaterial *mat = new TGeoMaterial("Al", 26.98,13,2.7);
28  TGeoMedium *med = new TGeoMedium("MED",1,mat);
29  TGeoVolume *top = gGeoManager->MakeBox("TOP",med,100,100,100);
30  gGeoManager->SetTopVolume(top);
31  TGeoVolume *vol = gGeoManager->MakeTorus("TORUS",med, 40,20,25,0,270);
32  top->AddNode(vol,1);
33  gGeoManager->CloseGeometry();
34  gGeoManager->SetNsegments(30);
35  top->Draw();
36  TView *view = gPad->GetView();
37  view->ShowAxis();
38 }
39 End_Macro
40 */
41 
42 #include "Riostream.h"
43 
44 #include "TGeoManager.h"
45 #include "TGeoVolume.h"
46 #include "TGeoTube.h"
47 #include "TVirtualGeoPainter.h"
48 #include "TGeoTorus.h"
49 #include "TVirtualPad.h"
50 #include "TBuffer3D.h"
51 #include "TBuffer3DTypes.h"
52 #include "TMath.h"
53 
54 ClassImp(TGeoTorus);
55 
56 ////////////////////////////////////////////////////////////////////////////////
57 /// Default constructor
58 
59 TGeoTorus::TGeoTorus()
60 {
61  SetShapeBit(TGeoShape::kGeoTorus);
62  fR = 0.0;
63  fRmin = 0.0;
64  fRmax = 0.0;
65  fPhi1 = 0.0;
66  fDphi = 0.0;
67 }
68 
69 ////////////////////////////////////////////////////////////////////////////////
70 /// Constructor without name.
71 
72 TGeoTorus::TGeoTorus(Double_t r, Double_t rmin, Double_t rmax, Double_t phi1, Double_t dphi)
73  :TGeoBBox(0, 0, 0)
74 {
75  SetShapeBit(TGeoShape::kGeoTorus);
76  SetTorusDimensions(r, rmin, rmax, phi1, dphi);
77  if ((fRmin<0) || (fRmax<0))
78  SetShapeBit(kGeoRunTimeShape);
79  ComputeBBox();
80 }
81 
82 ////////////////////////////////////////////////////////////////////////////////
83 /// Constructor with name.
84 
85 TGeoTorus::TGeoTorus(const char *name, Double_t r, Double_t rmin, Double_t rmax, Double_t phi1, Double_t dphi)
86  :TGeoBBox(name, 0, 0, 0)
87 {
88  SetShapeBit(TGeoShape::kGeoTorus);
89  SetTorusDimensions(r, rmin, rmax, phi1, dphi);
90  if ((fRmin<0) || (fRmax<0))
91  SetShapeBit(kGeoRunTimeShape);
92  ComputeBBox();
93 }
94 
95 ////////////////////////////////////////////////////////////////////////////////
96 /// Constructor based on an array of parameters.
97 /// - param[0] = R
98 /// - param[1] = Rmin
99 /// - param[2] = Rmax
100 /// - param[3] = Phi1
101 /// - param[4] = Dphi
102 
103 TGeoTorus::TGeoTorus(Double_t *param)
104  :TGeoBBox(0, 0, 0)
105 {
106  SetShapeBit(TGeoShape::kGeoTorus);
107  SetDimensions(param);
108  if (fRmin<0 || fRmax<0) SetShapeBit(kGeoRunTimeShape);
109  ComputeBBox();
110 }
111 
112 ////////////////////////////////////////////////////////////////////////////////
113 /// Computes capacity of the shape in [length^3]
114 
115 Double_t TGeoTorus::Capacity() const
116 {
117  Double_t capacity = (fDphi/180.)*TMath::Pi()*TMath::Pi()*fR*(fRmax*fRmax-fRmin*fRmin);
118  return capacity;
119 }
120 
121 ////////////////////////////////////////////////////////////////////////////////
122 /// Compute bounding box of the torus.
123 
124 void TGeoTorus::ComputeBBox()
125 {
126  fDZ = fRmax;
127  if (TGeoShape::IsSameWithinTolerance(fDphi,360)) {
128  fDX = fDY = fR+fRmax;
129  return;
130  }
131  Double_t xc[4];
132  Double_t yc[4];
133  xc[0] = (fR+fRmax)*TMath::Cos(fPhi1*TMath::DegToRad());
134  yc[0] = (fR+fRmax)*TMath::Sin(fPhi1*TMath::DegToRad());
135  xc[1] = (fR+fRmax)*TMath::Cos((fPhi1+fDphi)*TMath::DegToRad());
136  yc[1] = (fR+fRmax)*TMath::Sin((fPhi1+fDphi)*TMath::DegToRad());
137  xc[2] = (fR-fRmax)*TMath::Cos(fPhi1*TMath::DegToRad());
138  yc[2] = (fR-fRmax)*TMath::Sin(fPhi1*TMath::DegToRad());
139  xc[3] = (fR-fRmax)*TMath::Cos((fPhi1+fDphi)*TMath::DegToRad());
140  yc[3] = (fR-fRmax)*TMath::Sin((fPhi1+fDphi)*TMath::DegToRad());
141 
142  Double_t xmin = xc[TMath::LocMin(4, &xc[0])];
143  Double_t xmax = xc[TMath::LocMax(4, &xc[0])];
144  Double_t ymin = yc[TMath::LocMin(4, &yc[0])];
145  Double_t ymax = yc[TMath::LocMax(4, &yc[0])];
146  Double_t ddp = -fPhi1;
147  if (ddp<0) ddp+= 360;
148  if (ddp<=fDphi) xmax = fR+fRmax;
149  ddp = 90-fPhi1;
150  if (ddp<0) ddp+= 360;
151  if (ddp>360) ddp-=360;
152  if (ddp<=fDphi) ymax = fR+fRmax;
153  ddp = 180-fPhi1;
154  if (ddp<0) ddp+= 360;
155  if (ddp>360) ddp-=360;
156  if (ddp<=fDphi) xmin = -(fR+fRmax);
157  ddp = 270-fPhi1;
158  if (ddp<0) ddp+= 360;
159  if (ddp>360) ddp-=360;
160  if (ddp<=fDphi) ymin = -(fR+fRmax);
161  fOrigin[0] = (xmax+xmin)/2;
162  fOrigin[1] = (ymax+ymin)/2;
163  fOrigin[2] = 0;
164  fDX = (xmax-xmin)/2;
165  fDY = (ymax-ymin)/2;
166 }
167 
168 ////////////////////////////////////////////////////////////////////////////////
169 /// Compute normal to closest surface from POINT.
170 
171 void TGeoTorus::ComputeNormal(const Double_t *point, const Double_t *dir, Double_t *norm)
172 {
173  Double_t phi = TMath::ATan2(point[1],point[0]);
174  if (fDphi<360) {
175  Double_t phi1 = fPhi1*TMath::DegToRad();
176  Double_t phi2 = (fPhi1+fDphi)*TMath::DegToRad();
177  Double_t c1 = TMath::Cos(phi1);
178  Double_t s1 = TMath::Sin(phi1);
179  Double_t c2 = TMath::Cos(phi2);
180  Double_t s2 = TMath::Sin(phi2);
181 
182  Double_t daxis = Daxis(point,dir,0);
183  if ((fRmax-daxis)>1E-5) {
184  if (TGeoShape::IsSameWithinTolerance(fRmin,0) || (daxis-fRmin)>1E-5) {
185  TGeoShape::NormalPhi(point,dir,norm,c1,s1,c2,s2);
186  return;
187  }
188  }
189  }
190  Double_t r0[3];
191  r0[0] = fR*TMath::Cos(phi);
192  r0[1] = fR*TMath::Sin(phi);
193  r0[2] = 0;
194  Double_t normsq = 0;
195  for (Int_t i=0; i<3; i++) {
196  norm[i] = point[i] - r0[i];
197  normsq += norm[i]*norm[i];
198  }
199 
200  normsq = TMath::Sqrt(normsq);
201  norm[0] /= normsq;
202  norm[1] /= normsq;
203  norm[2] /= normsq;
204  if (dir[0]*norm[0]+dir[1]*norm[1]+dir[2]*norm[2] < 0) {
205  norm[0] = -norm[0];
206  norm[1] = -norm[1];
207  norm[2] = -norm[2];
208  }
209 }
210 
211 ////////////////////////////////////////////////////////////////////////////////
212 /// Test if point is inside the torus.
213 /// check phi range
214 
215 Bool_t TGeoTorus::Contains(const Double_t *point) const
216 {
217  if (!TGeoShape::IsSameWithinTolerance(fDphi,360)) {
218  Double_t phi = TMath::ATan2(point[1], point[0]) * TMath::RadToDeg();
219  if (phi < 0) phi+=360.0;
220  Double_t ddp = phi-fPhi1;
221  if (ddp<0) ddp+=360.;
222  if (ddp>fDphi) return kFALSE;
223  }
224  //check radius
225  Double_t rxy = TMath::Sqrt(point[0]*point[0]+point[1]*point[1]);
226  Double_t radsq = (rxy-fR)*(rxy-fR) + point[2]*point[2];
227  if (radsq<fRmin*fRmin) return kFALSE;
228  if (radsq>fRmax*fRmax) return kFALSE;
229  return kTRUE;
230 }
231 
232 ////////////////////////////////////////////////////////////////////////////////
233 /// Compute closest distance from point px,py to each vertex.
234 
235 Int_t TGeoTorus::DistancetoPrimitive(Int_t px, Int_t py)
236 {
237  Int_t n = gGeoManager->GetNsegments()+1;
238  Int_t numPoints = n*(n-1);
239  if (fRmin>0) numPoints *= 2;
240  else if (fDphi<360) numPoints += 2;
241  return ShapeDistancetoPrimitive(numPoints, px, py);
242 }
243 
244 ////////////////////////////////////////////////////////////////////////////////
245 /// Computes distance to axis of the torus from point pt + t*dir;
246 
247 Double_t TGeoTorus::Daxis(const Double_t *pt, const Double_t *dir, Double_t t) const
248 {
249  Double_t p[3];
250  for (Int_t i=0; i<3; i++) p[i] = pt[i]+t*dir[i];
251  Double_t rxy = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]);
252  return TMath::Sqrt((rxy-fR)*(rxy-fR)+p[2]*p[2]);
253 }
254 
255 ////////////////////////////////////////////////////////////////////////////////
256 /// Computes derivative w.r.t. t of the distance to axis of the torus from point pt + t*dir;
257 
258 Double_t TGeoTorus::DDaxis(const Double_t *pt, const Double_t *dir, Double_t t) const
259 {
260  Double_t p[3];
261  for (Int_t i=0; i<3; i++) p[i] = pt[i]+t*dir[i];
262  Double_t rxy = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]);
263  if (rxy<1E-4) return ((p[2]*dir[2]-fR*TMath::Sqrt(dir[0]*dir[0]+dir[1]*dir[1]))/TMath::Sqrt(fR*fR+p[2]*p[2]));
264  Double_t d = TMath::Sqrt((rxy-fR)*(rxy-fR)+p[2]*p[2]);
265  if (TGeoShape::IsSameWithinTolerance(d,0)) return 0.;
266  Double_t dd = (p[0]*dir[0]+p[1]*dir[1]+p[2]*dir[2] - (p[0]*dir[0]+p[1]*dir[1])*fR/rxy)/d;
267  return dd;
268 }
269 
270 ////////////////////////////////////////////////////////////////////////////////
271 /// Second derivative of distance to torus axis w.r.t t.
272 
273 Double_t TGeoTorus::DDDaxis(const Double_t *pt, const Double_t *dir, Double_t t) const
274 {
275  Double_t p[3];
276  for (Int_t i=0; i<3; i++) p[i] = pt[i]+t*dir[i];
277  Double_t rxy = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]);
278  if (rxy<1E-6) return 0;
279  Double_t daxis = TMath::Sqrt((rxy-fR)*(rxy-fR)+p[2]*p[2]);
280  if (TGeoShape::IsSameWithinTolerance(daxis,0)) return 0;
281  Double_t ddaxis = (p[0]*dir[0]+p[1]*dir[1]+p[2]*dir[2] - (p[0]*dir[0]+p[1]*dir[1])*fR/rxy)/daxis;
282  Double_t dddaxis = 1 - ddaxis*ddaxis - (1-dir[2]*dir[2])*fR/rxy +
283  fR*(p[0]*dir[0]+p[1]*dir[1])*(p[0]*dir[0]+p[1]*dir[1])/(rxy*rxy*rxy);
284  dddaxis /= daxis;
285  return dddaxis;
286 }
287 
288 ////////////////////////////////////////////////////////////////////////////////
289 /// Compute distance from inside point to surface of the torus.
290 
291 Double_t TGeoTorus::DistFromInside(const Double_t *point, const Double_t *dir, Int_t iact, Double_t step, Double_t *safe) const
292 {
293  if (iact<3 && safe) {
294  *safe = Safety(point, kTRUE);
295  if (iact==0) return TGeoShape::Big();
296  if ((iact==1) && (step<=*safe)) return TGeoShape::Big();
297  }
298  Double_t snext = TGeoShape::Big();
299  Bool_t hasphi = (fDphi<360)?kTRUE:kFALSE;
300  Bool_t hasrmin = (fRmin>0)?kTRUE:kFALSE;
301  Double_t dout = ToBoundary(point,dir,fRmax,kTRUE);
302 // Double_t dax = Daxis(point,dir,dout);
303  Double_t din = (hasrmin)?ToBoundary(point,dir,fRmin,kTRUE):TGeoShape::Big();
304  snext = TMath::Min(dout,din);
305  if (snext>1E10) return TGeoShape::Tolerance();
306  Double_t dphi = TGeoShape::Big();
307  if (hasphi) {
308  // Torus segment case.
309  Double_t c1,s1,c2,s2,cm,sm,cdfi;
310  Double_t phi1=fPhi1*TMath::DegToRad();
311  Double_t phi2=(fPhi1+fDphi)*TMath::DegToRad();
312  c1=TMath::Cos(phi1);
313  s1=TMath::Sin(phi1);
314  c2=TMath::Cos(phi2);
315  s2=TMath::Sin(phi2);
316  Double_t fio=0.5*(phi1+phi2);
317  cm=TMath::Cos(fio);
318  sm=TMath::Sin(fio);
319  cdfi = TMath::Cos(0.5*(phi2-phi1));
320  dphi = TGeoTubeSeg::DistFromInsideS(point,dir,fR-fRmax,fR+fRmax, fRmax, c1,s1,c2,s2,cm,sm,cdfi);
321  Double_t daxis = Daxis(point,dir,dphi);
322  if (daxis>=fRmin+1.E-8 && daxis<=fRmax-1.E-8) snext=TMath::Min(snext,dphi);
323  }
324  return snext;
325 }
326 
327 ////////////////////////////////////////////////////////////////////////////////
328 /// Compute distance from outside point to surface of the torus.
329 
330 Double_t TGeoTorus::DistFromOutside(const Double_t *point, const Double_t *dir, Int_t iact, Double_t step, Double_t *safe) const
331 {
332  if (iact<3 && safe) {
333  *safe = Safety(point, kFALSE);
334  if (iact==0) return TGeoShape::Big();
335  if ((iact==1) && (step<=*safe)) return TGeoShape::Big();
336  }
337 // Check if the bounding box is crossed within the requested distance
338  Double_t sdist = TGeoBBox::DistFromOutside(point,dir, fDX, fDY, fDZ, fOrigin, step);
339  if (sdist>=step) return TGeoShape::Big();
340  Double_t daxis;
341  Bool_t hasphi = (fDphi<360)?kTRUE:kFALSE;
342 // Bool_t hasrmin = (fRmin>0)?kTRUE:kFALSE;
343  Double_t c1=0,s1=0,c2=0,s2=0,cm=0,sm=0,cdfi=0;
344  Bool_t inphi = kFALSE;
345  Double_t phi, ddp, phi1,phi2,fio;
346  Double_t rxy2,dd;
347  Double_t snext;
348  Double_t pt[3];
349  Int_t i;
350 
351  if (hasphi) {
352  // Torus segment case.
353  phi=TMath::ATan2(point[1], point[0])*TMath::RadToDeg();;
354  if (phi<0) phi+=360;
355  ddp = phi-fPhi1;
356  if (ddp<0) ddp+=360;;
357  if (ddp<=fDphi) inphi=kTRUE;
358  phi1=fPhi1*TMath::DegToRad();
359  phi2=(fPhi1+fDphi)*TMath::DegToRad();
360  c1=TMath::Cos(phi1);
361  s1=TMath::Sin(phi1);
362  c2=TMath::Cos(phi2);
363  s2=TMath::Sin(phi2);
364  fio=0.5*(phi1+phi2);
365  cm=TMath::Cos(fio);
366  sm=TMath::Sin(fio);
367  cdfi=TMath::Cos(0.5*(phi2-phi1));
368  }
369  // Check if we are inside or outside the bounding ring.
370  Bool_t inbring = kFALSE;
371  if (TMath::Abs(point[2]) <= fRmax) {
372  rxy2 = point[0]*point[0]+point[1]*point[1];
373  if ((rxy2>=(fR-fRmax)*(fR-fRmax)) && (rxy2<=(fR+fRmax)*(fR+fRmax))) {
374  if (!hasphi || inphi) inbring=kTRUE;
375  }
376  }
377 
378  // If outside the ring, compute distance to it.
379  Double_t dring = TGeoShape::Big();
380  Double_t eps = 1.E-8;
381  snext = 0;
382  daxis = -1;
383  memcpy(pt,point,3*sizeof(Double_t));
384  if (!inbring) {
385  if (hasphi) dring = TGeoTubeSeg::DistFromOutsideS(point,dir,TMath::Max(0.,fR-fRmax-eps),fR+fRmax+eps, fRmax+eps, c1,s1,c2,s2,cm,sm,cdfi);
386  else dring = TGeoTube::DistFromOutsideS(point,dir,TMath::Max(0.,fR-fRmax-eps),fR+fRmax+eps, fRmax+eps);
387  // If not crossing it, return BIG.
388  if (dring>1E10) return TGeoShape::Big();
389  snext = dring;
390  // Check if the crossing is due to phi.
391  daxis = Daxis(point,dir,snext);
392  if (daxis>=fRmin && daxis<fRmax) return snext;
393  // Not a phi crossing -> propagate until we cross the ring.
394  for (i=0; i<3; i++) pt[i] = point[i]+snext*dir[i];
395  }
396  // Point pt is inside the bounding ring, no phi crossing so far.
397  // Check if we are in the hole.
398  if (daxis<0) daxis = Daxis(pt,dir,0);
399  if (daxis<fRmin+1.E-8) {
400  // We are in the hole. Check if we came from outside.
401  if (snext>0) {
402  // we can cross either the inner torus or exit the other hole.
403  snext += 0.1*eps;
404  for (i=0; i<3; i++) pt[i] += 0.1*eps*dir[i];
405  }
406  // We are in the hole from the beginning.
407  // find first crossing with inner torus
408  dd = ToBoundary(pt,dir, fRmin,kFALSE);
409  // find exit distance from inner bounding ring
410  if (hasphi) dring = TGeoTubeSeg::DistFromInsideS(pt,dir,fR-fRmin,fR+fRmin, fRmin, c1,s1,c2,s2,cm,sm,cdfi);
411  else dring = TGeoTube::DistFromInsideS(pt,dir,fR-fRmin,fR+fRmin, fRmin);
412  if (dd<dring) return (snext+dd);
413  // we were exiting a hole inside phi hole
414  snext += dring+ eps;
415  for (i=0; i<3; i++) pt[i] = point[i] + snext*dir[i];
416  snext += DistFromOutside(pt,dir,3);
417  return snext;
418  }
419  // We are inside the outer ring, having daxis>fRmax
420  // Compute distance to exit the bounding ring (again)
421  if (snext>0) {
422  // we can cross either the inner torus or exit the other hole.
423  snext += 0.1*eps;
424  for (i=0; i<3; i++) pt[i] += 0.1*eps*dir[i];
425  }
426  // Check intersection with outer torus
427  dd = ToBoundary(pt, dir, fRmax, kFALSE);
428  if (hasphi) dring = TGeoTubeSeg::DistFromInsideS(pt,dir,TMath::Max(0.,fR-fRmax-eps),fR+fRmax+eps, fRmax+eps, c1,s1,c2,s2,cm,sm,cdfi);
429  else dring = TGeoTube::DistFromInsideS(pt,dir,TMath::Max(0.,fR-fRmax-eps),fR+fRmax+eps, fRmax+eps);
430  if (dd<dring) {
431  snext += dd;
432  return snext;
433  }
434  // We are exiting the bounding ring before crossing the torus -> propagate
435  snext += dring+eps;
436  for (i=0; i<3; i++) pt[i] = point[i] + snext*dir[i];
437  snext += DistFromOutside(pt,dir,3);
438  return snext;
439 }
440 
441 ////////////////////////////////////////////////////////////////////////////////
442 /// Divide this torus shape belonging to volume "voldiv" into ndiv volumes
443 /// called divname, from start position with the given step.
444 
445 TGeoVolume *TGeoTorus::Divide(TGeoVolume * /*voldiv*/, const char * /*divname*/, Int_t /*iaxis*/, Int_t /*ndiv*/,
446  Double_t /*start*/, Double_t /*step*/)
447 {
448  return 0;
449 }
450 
451 ////////////////////////////////////////////////////////////////////////////////
452 /// Returns name of axis IAXIS.
453 
454 const char *TGeoTorus::GetAxisName(Int_t iaxis) const
455 {
456  switch (iaxis) {
457  case 1:
458  return "R";
459  case 2:
460  return "PHI";
461  case 3:
462  return "Z";
463  default:
464  return "UNDEFINED";
465  }
466 }
467 
468 ////////////////////////////////////////////////////////////////////////////////
469 /// Get range of shape for a given axis.
470 
471 Double_t TGeoTorus::GetAxisRange(Int_t iaxis, Double_t &xlo, Double_t &xhi) const
472 {
473  xlo = 0;
474  xhi = 0;
475  Double_t dx = 0;
476  switch (iaxis) {
477  case 1:
478  xlo = fRmin;
479  xhi = fRmax;
480  dx = xhi-xlo;
481  return dx;
482  case 2:
483  xlo = fPhi1;
484  xhi = fPhi1+fDphi;
485  dx = fDphi;
486  return dx;
487  case 3:
488  dx = 0;
489  return dx;
490  }
491  return dx;
492 }
493 
494 ////////////////////////////////////////////////////////////////////////////////
495 /// Fill vector param[4] with the bounding cylinder parameters. The order
496 /// is the following : Rmin, Rmax, Phi1, Phi2, dZ
497 
498 void TGeoTorus::GetBoundingCylinder(Double_t *param) const
499 {
500  param[0] = (fR-fRmax); // Rmin
501  param[1] = (fR+fRmax); // Rmax
502  param[2] = fPhi1; // Phi1
503  param[3] = fPhi1+fDphi; // Phi2
504 }
505 
506 ////////////////////////////////////////////////////////////////////////////////
507 /// Create a shape fitting the mother.
508 
509 TGeoShape *TGeoTorus::GetMakeRuntimeShape(TGeoShape * /*mother*/, TGeoMatrix * /*mat*/) const
510 {
511  if (!TestShapeBit(kGeoRunTimeShape)) return 0;
512  Error("GetMakeRuntimeShape", "parametrized toruses not supported");
513  return 0;
514 }
515 
516 ////////////////////////////////////////////////////////////////////////////////
517 /// print shape parameters
518 
519 void TGeoTorus::InspectShape() const
520 {
521  printf("*** Shape %s: TGeoTorus ***\n", GetName());
522  printf(" R = %11.5f\n", fR);
523  printf(" Rmin = %11.5f\n", fRmin);
524  printf(" Rmax = %11.5f\n", fRmax);
525  printf(" Phi1 = %11.5f\n", fPhi1);
526  printf(" Dphi = %11.5f\n", fDphi);
527  printf(" Bounding box:\n");
528  TGeoBBox::InspectShape();
529 }
530 
531 ////////////////////////////////////////////////////////////////////////////////
532 /// Creates a TBuffer3D describing *this* shape.
533 /// Coordinates are in local reference frame.
534 
535 TBuffer3D *TGeoTorus::MakeBuffer3D() const
536 {
537  Int_t n = gGeoManager->GetNsegments()+1;
538  Int_t nbPnts = n*(n-1);
539  Bool_t hasrmin = (GetRmin()>0)?kTRUE:kFALSE;
540  Bool_t hasphi = (GetDphi()<360)?kTRUE:kFALSE;
541  if (hasrmin) nbPnts *= 2;
542  else if (hasphi) nbPnts += 2;
543 
544  Int_t nbSegs = (2*n-1)*(n-1);
545  Int_t nbPols = (n-1)*(n-1);
546  if (hasrmin) {
547  nbSegs += (2*n-1)*(n-1);
548  nbPols += (n-1)*(n-1);
549  }
550  if (hasphi) {
551  nbSegs += 2*(n-1);
552  nbPols += 2*(n-1);
553  }
554 
555  TBuffer3D* buff = new TBuffer3D(TBuffer3DTypes::kGeneric,
556  nbPnts, 3*nbPnts, nbSegs, 3*nbSegs, nbPols, 6*nbPols);
557  if (buff)
558  {
559  SetPoints(buff->fPnts);
560  SetSegsAndPols(*buff);
561  }
562 
563  return buff;
564 }
565 
566 ////////////////////////////////////////////////////////////////////////////////
567 /// Fill TBuffer3D structure for segments and polygons.
568 
569 void TGeoTorus::SetSegsAndPols(TBuffer3D &buff) const
570 {
571  Int_t i, j;
572  Int_t n = gGeoManager->GetNsegments()+1;
573  Int_t nbPnts = n*(n-1);
574  Int_t indx, indp, startcap=0;
575  Bool_t hasrmin = (GetRmin()>0)?kTRUE:kFALSE;
576  Bool_t hasphi = (GetDphi()<360)?kTRUE:kFALSE;
577  if (hasrmin) nbPnts *= 2;
578  else if (hasphi) nbPnts += 2;
579  Int_t c = GetBasicColor();
580 
581  indp = n*(n-1); // start index for points on inner surface
582  memset(buff.fSegs, 0, buff.NbSegs()*3*sizeof(Int_t));
583 
584  // outer surface phi circles = n*(n-1) -> [0, n*(n-1) -1]
585  // connect point j with point j+1 on same row
586  indx = 0;
587  for (i = 0; i < n; i++) { // rows [0,n-1]
588  for (j = 0; j < n-1; j++) { // points on a row [0, n-2]
589  buff.fSegs[indx+(i*(n-1)+j)*3] = c;
590  buff.fSegs[indx+(i*(n-1)+j)*3+1] = i*(n-1)+j; // j on row i
591  buff.fSegs[indx+(i*(n-1)+j)*3+2] = i*(n-1)+((j+1)%(n-1)); // j+1 on row i
592  }
593  }
594  indx += 3*n*(n-1);
595  // outer surface generators = (n-1)*(n-1) -> [n*(n-1), (2*n-1)*(n-1) -1]
596  // connect point j on row i with point j on row i+1
597  for (i = 0; i < n-1; i++) { // rows [0, n-2]
598  for (j = 0; j < n-1; j++) { // points on a row [0, n-2]
599  buff.fSegs[indx+(i*(n-1)+j)*3] = c;
600  buff.fSegs[indx+(i*(n-1)+j)*3+1] = i*(n-1)+j; // j on row i
601  buff.fSegs[indx+(i*(n-1)+j)*3+2] = (i+1)*(n-1)+j; // j on row i+1
602  }
603  }
604  indx += 3*(n-1)*(n-1);
605  startcap = (2*n-1)*(n-1);
606 
607  if (hasrmin) {
608  // inner surface phi circles = n*(n-1) -> [(2*n-1)*(n-1), (3*n-1)*(n-1) -1]
609  // connect point j with point j+1 on same row
610  for (i = 0; i < n; i++) { // rows [0, n-1]
611  for (j = 0; j < n-1; j++) { // points on a row [0, n-2]
612  buff.fSegs[indx+(i*(n-1)+j)*3] = c; // lighter color
613  buff.fSegs[indx+(i*(n-1)+j)*3+1] = indp + i*(n-1)+j; // j on row i
614  buff.fSegs[indx+(i*(n-1)+j)*3+2] = indp + i*(n-1)+((j+1)%(n-1)); // j+1 on row i
615  }
616  }
617  indx += 3*n*(n-1);
618  // inner surface generators = (n-1)*n -> [(3*n-1)*(n-1), (4*n-2)*(n-1) -1]
619  // connect point j on row i with point j on row i+1
620  for (i = 0; i < n-1; i++) { // rows [0, n-2]
621  for (j = 0; j < n-1; j++) { // points on a row [0, n-2]
622  buff.fSegs[indx+(i*(n-1)+j)*3] = c; // lighter color
623  buff.fSegs[indx+(i*(n-1)+j)*3+1] = indp + i*(n-1)+j; // j on row i
624  buff.fSegs[indx+(i*(n-1)+j)*3+2] = indp + (i+1)*(n-1)+j; // j on row i+1
625  }
626  }
627  indx += 3*(n-1)*(n-1);
628  startcap = (4*n-2)*(n-1);
629  }
630 
631  if (hasphi) {
632  if (hasrmin) {
633  // endcaps = 2*(n-1) -> [(4*n-2)*(n-1), 4*n*(n-1)-1]
634  i = 0;
635  for (j = 0; j < n-1; j++) {
636  buff.fSegs[indx+j*3] = c+1;
637  buff.fSegs[indx+j*3+1] = (n-1)*i+j; // outer j on row 0
638  buff.fSegs[indx+j*3+2] = indp+(n-1)*i+j; // inner j on row 0
639  }
640  indx += 3*(n-1);
641  i = n-1;
642  for (j = 0; j < n-1; j++) {
643  buff.fSegs[indx+j*3] = c+1;
644  buff.fSegs[indx+j*3+1] = (n-1)*i+j; // outer j on row n-1
645  buff.fSegs[indx+j*3+2] = indp+(n-1)*i+j; // inner j on row n-1
646  }
647  indx += 3*(n-1);
648  } else {
649  i = 0;
650  for (j = 0; j < n-1; j++) {
651  buff.fSegs[indx+j*3] = c+1;
652  buff.fSegs[indx+j*3+1] = (n-1)*i+j; // outer j on row 0
653  buff.fSegs[indx+j*3+2] = n*(n-1); // center of first endcap
654  }
655  indx += 3*(n-1);
656  i = n-1;
657  for (j = 0; j < n-1; j++) {
658  buff.fSegs[indx+j*3] = c+1;
659  buff.fSegs[indx+j*3+1] = (n-1)*i+j; // outer j on row n-1
660  buff.fSegs[indx+j*3+2] = n*(n-1)+1; // center of second endcap
661  }
662  indx += 3*(n-1);
663  }
664  }
665 
666  indx = 0;
667  memset(buff.fPols, 0, buff.NbPols()*6*sizeof(Int_t));
668 
669  // outer surface = (n-1)*(n-1) -> [0, (n-1)*(n-1)-1]
670  // normal pointing out
671  for (i=0; i<n-1; i++) {
672  for (j=0; j<n-1; j++) {
673  buff.fPols[indx++] = c;
674  buff.fPols[indx++] = 4;
675  buff.fPols[indx++] = n*(n-1)+(n-1)*i+((j+1)%(n-1)); // generator j+1 on outer row i
676  buff.fPols[indx++] = (n-1)*(i+1)+j; // seg j on outer row i+1
677  buff.fPols[indx++] = n*(n-1)+(n-1)*i+j; // generator j on outer row i
678  buff.fPols[indx++] = (n-1)*i+j; // seg j on outer row i
679  }
680  }
681  if (hasrmin) {
682  indp = (2*n-1)*(n-1); // start index of inner segments
683  // inner surface = (n-1)*(n-1) -> [(n-1)*(n-1), 2*(n-1)*(n-1)-1]
684  // normal pointing out
685  for (i=0; i<n-1; i++) {
686  for (j=0; j<n-1; j++) {
687  buff.fPols[indx++] = c;
688  buff.fPols[indx++] = 4;
689  buff.fPols[indx++] = indp+n*(n-1)+(n-1)*i+j; // generator j on inner row i
690  buff.fPols[indx++] = indp+(n-1)*(i+1)+j; // seg j on inner row i+1
691  buff.fPols[indx++] = indp+n*(n-1)+(n-1)*i+((j+1)%(n-1)); // generator j+1 on inner r>
692  buff.fPols[indx++] = indp+(n-1)*i+j; // seg j on inner row i
693  }
694  }
695  }
696  if (hasphi) {
697  // endcaps = 2*(n-1) -> [2*(n-1)*(n-1), 2*n*(n-1)-1]
698  i=0; // row 0
699  Int_t np = (hasrmin)?4:3;
700  for (j=0; j<n-1; j++) {
701  buff.fPols[indx++] = c+1;
702  buff.fPols[indx++] = np;
703  buff.fPols[indx++] = j; // seg j on outer row 0 a
704  buff.fPols[indx++] = startcap+j; // endcap j on row 0 d
705  if(hasrmin)
706  buff.fPols[indx++] = indp+j; // seg j on inner row 0 c
707  buff.fPols[indx++] = startcap+((j+1)%(n-1)); // endcap j+1 on row 0 b
708  }
709 
710  i=n-1; // row n-1
711  for (j=0; j<n-1; j++) {
712  buff.fPols[indx++] = c+1;
713  buff.fPols[indx++] = np;
714  buff.fPols[indx++] = (n-1)*i+j; // seg j on outer row n-1 a
715  buff.fPols[indx++] = startcap+(n-1)+((j+1)%(n-1)); // endcap j+1 on row n-1 d
716  if (hasrmin)
717  buff.fPols[indx++] = indp+(n-1)*i+j; // seg j on inner row n-1 c
718  buff.fPols[indx++] = startcap+(n-1)+j; // endcap j on row n-1 b
719  }
720  }
721 }
722 
723 ////////////////////////////////////////////////////////////////////////////////
724 /// computes the closest distance from given point to this shape, according
725 /// to option. The matching point on the shape is stored in spoint.
726 
727 Double_t TGeoTorus::Safety(const Double_t *point, Bool_t in) const
728 {
729  Double_t saf[2];
730  Int_t i;
731  Double_t rxy = TMath::Sqrt(point[0]*point[0]+point[1]*point[1]);
732  Double_t rad = TMath::Sqrt((rxy-fR)*(rxy-fR) + point[2]*point[2]);
733  saf[0] = rad-fRmin;
734  saf[1] = fRmax-rad;
735  if (TGeoShape::IsSameWithinTolerance(fDphi,360)) {
736  if (in) return TMath::Min(saf[0],saf[1]);
737  for (i=0; i<2; i++) saf[i]=-saf[i];
738  return TMath::Max(saf[0], saf[1]);
739  }
740 
741  Double_t safphi = TGeoShape::SafetyPhi(point,in,fPhi1, fPhi1+fDphi);
742  Double_t safe = TGeoShape::Big();
743  if (in) {
744  safe = TMath::Min(saf[0], saf[1]);
745  return TMath::Min(safe, safphi);
746  }
747  for (i=0; i<2; i++) saf[i]=-saf[i];
748  safe = TMath::Max(saf[0], saf[1]);
749  return TMath::Max(safe, safphi);
750 }
751 
752 ////////////////////////////////////////////////////////////////////////////////
753 /// Save a primitive as a C++ statement(s) on output stream "out".
754 
755 void TGeoTorus::SavePrimitive(std::ostream &out, Option_t * /*option*/ /*= ""*/)
756 {
757  if (TObject::TestBit(kGeoSavePrimitive)) return;
758  out << " // Shape: " << GetName() << " type: " << ClassName() << std::endl;
759  out << " r = " << fR << ";" << std::endl;
760  out << " rmin = " << fRmin << ";" << std::endl;
761  out << " rmax = " << fRmax << ";" << std::endl;
762  out << " phi1 = " << fPhi1 << ";" << std::endl;
763  out << " dphi = " << fDphi << ";" << std::endl;
764  out << " TGeoShape *" << GetPointerName() << " = new TGeoTorus(\"" << GetName() << "\",r,rmin,rmax,phi1,dphi);" << std::endl;
765  TObject::SetBit(TGeoShape::kGeoSavePrimitive);
766 }
767 
768 ////////////////////////////////////////////////////////////////////////////////
769 /// Set torus dimensions.
770 
771 void TGeoTorus::SetTorusDimensions(Double_t r, Double_t rmin, Double_t rmax,
772  Double_t phi1, Double_t dphi)
773 {
774  fR = r;
775  fRmin = rmin;
776  fRmax = rmax;
777  fPhi1 = phi1;
778  if (fPhi1<0) fPhi1+=360.;
779  fDphi = dphi;
780 }
781 
782 ////////////////////////////////////////////////////////////////////////////////
783 /// Set torus dimensions starting from a list.
784 
785 void TGeoTorus::SetDimensions(Double_t *param)
786 {
787  SetTorusDimensions(param[0], param[1], param[2], param[3], param[4]);
788 }
789 
790 ////////////////////////////////////////////////////////////////////////////////
791 /// Create torus mesh points
792 
793 void TGeoTorus::SetPoints(Double_t *points) const
794 {
795  if (!points) return;
796  Int_t n = gGeoManager->GetNsegments()+1;
797  Double_t phin, phout;
798  Double_t dpin = 360./(n-1);
799  Double_t dpout = fDphi/(n-1);
800  Double_t co,so,ci,si;
801  Bool_t havermin = (fRmin<TGeoShape::Tolerance())?kFALSE:kTRUE;
802  Int_t i,j;
803  Int_t indx = 0;
804  // loop outer mesh -> n*(n-1) points [0, 3*n*(n-1)-1]
805  for (i=0; i<n; i++) {
806  phout = (fPhi1+i*dpout)*TMath::DegToRad();
807  co = TMath::Cos(phout);
808  so = TMath::Sin(phout);
809  for (j=0; j<n-1; j++) {
810  phin = j*dpin*TMath::DegToRad();
811  ci = TMath::Cos(phin);
812  si = TMath::Sin(phin);
813  points[indx++] = (fR+fRmax*ci)*co;
814  points[indx++] = (fR+fRmax*ci)*so;
815  points[indx++] = fRmax*si;
816  }
817  }
818 
819  if (havermin) {
820  // loop inner mesh -> n*(n-1) points [3*n*(n-1), 6*n*(n-1)]
821  for (i=0; i<n; i++) {
822  phout = (fPhi1+i*dpout)*TMath::DegToRad();
823  co = TMath::Cos(phout);
824  so = TMath::Sin(phout);
825  for (j=0; j<n-1; j++) {
826  phin = j*dpin*TMath::DegToRad();
827  ci = TMath::Cos(phin);
828  si = TMath::Sin(phin);
829  points[indx++] = (fR+fRmin*ci)*co;
830  points[indx++] = (fR+fRmin*ci)*so;
831  points[indx++] = fRmin*si;
832  }
833  }
834  } else {
835  if (fDphi<360.) {
836  // just add extra 2 points on the centers of the 2 phi cuts [3*n*n, 3*n*n+1]
837  points[indx++] = fR*TMath::Cos(fPhi1*TMath::DegToRad());
838  points[indx++] = fR*TMath::Sin(fPhi1*TMath::DegToRad());
839  points[indx++] = 0;
840  points[indx++] = fR*TMath::Cos((fPhi1+fDphi)*TMath::DegToRad());
841  points[indx++] = fR*TMath::Sin((fPhi1+fDphi)*TMath::DegToRad());
842  points[indx++] = 0;
843  }
844  }
845 }
846 
847 ////////////////////////////////////////////////////////////////////////////////
848 /// Create torus mesh points
849 
850 void TGeoTorus::SetPoints(Float_t *points) const
851 {
852  if (!points) return;
853  Int_t n = gGeoManager->GetNsegments()+1;
854  Double_t phin, phout;
855  Double_t dpin = 360./(n-1);
856  Double_t dpout = fDphi/(n-1);
857  Double_t co,so,ci,si;
858  Bool_t havermin = (fRmin<TGeoShape::Tolerance())?kFALSE:kTRUE;
859  Int_t i,j;
860  Int_t indx = 0;
861  // loop outer mesh -> n*(n-1) points [0, 3*n*(n-1)-1]
862  // plane i = 0, n-1 point j = 0, n-1 ipoint = n*i + j
863  for (i=0; i<n; i++) {
864  phout = (fPhi1+i*dpout)*TMath::DegToRad();
865  co = TMath::Cos(phout);
866  so = TMath::Sin(phout);
867  for (j=0; j<n-1; j++) {
868  phin = j*dpin*TMath::DegToRad();
869  ci = TMath::Cos(phin);
870  si = TMath::Sin(phin);
871  points[indx++] = (fR+fRmax*ci)*co;
872  points[indx++] = (fR+fRmax*ci)*so;
873  points[indx++] = fRmax*si;
874  }
875  }
876 
877  if (havermin) {
878  // loop inner mesh -> n*(n-1) points [3*n*(n-1), 6*n*(n-1)]
879  // plane i = 0, n-1 point j = 0, n-1 ipoint = n*n + n*i + j
880  for (i=0; i<n; i++) {
881  phout = (fPhi1+i*dpout)*TMath::DegToRad();
882  co = TMath::Cos(phout);
883  so = TMath::Sin(phout);
884  for (j=0; j<n-1; j++) {
885  phin = j*dpin*TMath::DegToRad();
886  ci = TMath::Cos(phin);
887  si = TMath::Sin(phin);
888  points[indx++] = (fR+fRmin*ci)*co;
889  points[indx++] = (fR+fRmin*ci)*so;
890  points[indx++] = fRmin*si;
891  }
892  }
893  } else {
894  if (fDphi<360.) {
895  // just add extra 2 points on the centers of the 2 phi cuts [n*n, n*n+1]
896  // ip1 = n*(n-1) + 0;
897  // ip2 = n*(n-1) + 1
898  points[indx++] = fR*TMath::Cos(fPhi1*TMath::DegToRad());
899  points[indx++] = fR*TMath::Sin(fPhi1*TMath::DegToRad());
900  points[indx++] = 0;
901  points[indx++] = fR*TMath::Cos((fPhi1+fDphi)*TMath::DegToRad());
902  points[indx++] = fR*TMath::Sin((fPhi1+fDphi)*TMath::DegToRad());
903  points[indx++] = 0;
904  }
905  }
906 }
907 
908 ////////////////////////////////////////////////////////////////////////////////
909 /// Return number of vertices of the mesh representation
910 
911 Int_t TGeoTorus::GetNmeshVertices() const
912 {
913  Int_t n = gGeoManager->GetNsegments()+1;
914  Int_t numPoints = n*(n-1);
915  if (fRmin>TGeoShape::Tolerance()) numPoints *= 2;
916  else if (fDphi<360.) numPoints += 2;
917  return numPoints;
918 }
919 
920 ////////////////////////////////////////////////////////////////////////////////
921 /// fill size of this 3-D object
922 
923 void TGeoTorus::Sizeof3D() const
924 {
925 }
926 
927 ////////////////////////////////////////////////////////////////////////////////
928 /// Find real solutions of the cubic equation : x^3 + a*x^2 + b*x + c = 0
929 /// Input: a,b,c
930 /// Output: x[3] real solutions
931 /// Returns number of real solutions (1 or 3)
932 
933 Int_t TGeoTorus::SolveCubic(Double_t a, Double_t b, Double_t c, Double_t *x) const
934 {
935  const Double_t ott = 1./3.;
936  const Double_t sq3 = TMath::Sqrt(3.);
937  Int_t ireal = 1;
938  Double_t p = b-a*a*ott;
939  Double_t q = c-a*b*ott+2.*a*a*a*ott*ott*ott;
940  Double_t delta = 4*p*p*p+27*q*q;
941 // Double_t y1r, y1i, y2r, y2i;
942  Double_t t,u;
943  if (delta>=0) {
944  delta = TMath::Sqrt(delta);
945  t = (-3*q*sq3+delta)/(6*sq3);
946  u = (3*q*sq3+delta)/(6*sq3);
947  x[0] = TMath::Sign(1.,t)*TMath::Power(TMath::Abs(t),ott)-
948  TMath::Sign(1.,u)*TMath::Power(TMath::Abs(u),ott)-a*ott;
949  } else {
950  delta = TMath::Sqrt(-delta);
951  t = -0.5*q;
952  u = delta/(6*sq3);
953  x[0] = 2.*TMath::Power(t*t+u*u,0.5*ott) * TMath::Cos(ott*TMath::ATan2(u,t));
954  x[0] -= a*ott;
955  }
956 
957  t = x[0]*x[0]+a*x[0]+b;
958  u = a+x[0];
959  delta = u*u-4.*t;
960  if (delta>=0) {
961  ireal = 3;
962  delta = TMath::Sqrt(delta);
963  x[1] = 0.5*(-u-delta);
964  x[2] = 0.5*(-u+delta);
965  }
966  return ireal;
967 }
968 
969 ////////////////////////////////////////////////////////////////////////////////
970 /// Find real solutions of the quartic equation : x^4 + a*x^3 + b*x^2 + c*x + d = 0
971 /// Input: a,b,c,d
972 /// Output: x[4] - real solutions
973 /// Returns number of real solutions (0 to 3)
974 
975 Int_t TGeoTorus::SolveQuartic(Double_t a, Double_t b, Double_t c, Double_t d, Double_t *x) const
976 {
977  Double_t e = b-3.*a*a/8.;
978  Double_t f = c+a*a*a/8.-0.5*a*b;
979  Double_t g = d-3.*a*a*a*a/256. + a*a*b/16. - a*c/4.;
980  Double_t xx[4];
981  Int_t ind[4];
982  Double_t delta;
983  Double_t h=0;
984  Int_t ireal = 0;
985  Int_t i;
986  if (TGeoShape::IsSameWithinTolerance(f,0)) {
987  delta = e*e-4.*g;
988  if (delta<0) return 0;
989  delta = TMath::Sqrt(delta);
990  h = 0.5*(-e-delta);
991  if (h>=0) {
992  h = TMath::Sqrt(h);
993  x[ireal++] = -h-0.25*a;
994  x[ireal++] = h-0.25*a;
995  }
996  h = 0.5*(-e+delta);
997  if (h>=0) {
998  h = TMath::Sqrt(h);
999  x[ireal++] = -h-0.25*a;
1000  x[ireal++] = h-0.25*a;
1001  }
1002  if (ireal>0) {
1003  TMath::Sort(ireal, x, ind,kFALSE);
1004  for (i=0; i<ireal; i++) xx[i] = x[ind[i]];
1005  memcpy(x,xx,ireal*sizeof(Double_t));
1006  }
1007  return ireal;
1008  }
1009 
1010  if (TGeoShape::IsSameWithinTolerance(g,0)) {
1011  x[ireal++] = -0.25*a;
1012  ind[0] = SolveCubic(0,e,f,xx);
1013  for (i=0; i<ind[0]; i++) x[ireal++] = xx[i]-0.25*a;
1014  if (ireal>0) {
1015  TMath::Sort(ireal, x, ind,kFALSE);
1016  for (i=0; i<ireal; i++) xx[i] = x[ind[i]];
1017  memcpy(x,xx,ireal*sizeof(Double_t));
1018  }
1019  return ireal;
1020  }
1021 
1022 
1023  ireal = SolveCubic(2.*e, e*e-4.*g, -f*f, xx);
1024  if (ireal==1) {
1025  if (xx[0]<=0) return 0;
1026  h = TMath::Sqrt(xx[0]);
1027  } else {
1028  // 3 real solutions of the cubic
1029  for (i=0; i<3; i++) {
1030  h = xx[i];
1031  if (h>=0) break;
1032  }
1033  if (h<=0) return 0;
1034  h = TMath::Sqrt(h);
1035  }
1036  Double_t j = 0.5*(e+h*h-f/h);
1037  ireal = 0;
1038  delta = h*h-4.*j;
1039  if (delta>=0) {
1040  delta = TMath::Sqrt(delta);
1041  x[ireal++] = 0.5*(-h-delta)-0.25*a;
1042  x[ireal++] = 0.5*(-h+delta)-0.25*a;
1043  }
1044  delta = h*h-4.*g/j;
1045  if (delta>=0) {
1046  delta = TMath::Sqrt(delta);
1047  x[ireal++] = 0.5*(h-delta)-0.25*a;
1048  x[ireal++] = 0.5*(h+delta)-0.25*a;
1049  }
1050  if (ireal>0) {
1051  TMath::Sort(ireal, x, ind,kFALSE);
1052  for (i=0; i<ireal; i++) xx[i] = x[ind[i]];
1053  memcpy(x,xx,ireal*sizeof(Double_t));
1054  }
1055  return ireal;
1056 }
1057 
1058 ////////////////////////////////////////////////////////////////////////////////
1059 /// Returns distance to the surface or the torus (fR,r) from a point, along
1060 /// a direction. Point is close enough to the boundary so that the distance
1061 /// to the torus is decreasing while moving along the given direction.
1062 
1063 Double_t TGeoTorus::ToBoundary(const Double_t *pt, const Double_t *dir, Double_t r, Bool_t in) const
1064 {
1065  // Compute coefficients of the quartic
1066  Double_t s = TGeoShape::Big();
1067  Double_t tol = TGeoShape::Tolerance();
1068  Double_t r0sq = pt[0]*pt[0]+pt[1]*pt[1]+pt[2]*pt[2];
1069  Double_t rdotn = pt[0]*dir[0]+pt[1]*dir[1]+pt[2]*dir[2];
1070  Double_t rsumsq = fR*fR+r*r;
1071  Double_t a = 4.*rdotn;
1072  Double_t b = 2.*(r0sq+2.*rdotn*rdotn-rsumsq+2.*fR*fR*dir[2]*dir[2]);
1073  Double_t c = 4.*(r0sq*rdotn-rsumsq*rdotn+2.*fR*fR*pt[2]*dir[2]);
1074  Double_t d = r0sq*r0sq-2.*r0sq*rsumsq+4.*fR*fR*pt[2]*pt[2]+(fR*fR-r*r)*(fR*fR-r*r);
1075 
1076  Double_t x[4],y[4];
1077  Int_t nsol = 0;
1078 
1079  if (TMath::Abs(dir[2])<1E-3 && TMath::Abs(pt[2])<0.1*r) {
1080  Double_t r0 = fR - TMath::Sqrt((r-pt[2])*(r+pt[2]));
1081  Double_t b0 = (pt[0]*dir[0]+pt[1]*dir[1])/(dir[0]*dir[0]+dir[1]*dir[1]);
1082  Double_t c0 = (pt[0]*pt[0] + (pt[1]-r0)*(pt[1]+r0))/(dir[0]*dir[0]+dir[1]*dir[1]);
1083  Double_t delta = b0*b0-c0;
1084  if (delta>0) {
1085  y[nsol] = -b0-TMath::Sqrt(delta);
1086  if (y[nsol]>-tol) nsol++;
1087  y[nsol] = -b0+TMath::Sqrt(delta);
1088  if (y[nsol]>-tol) nsol++;
1089  }
1090  r0 = fR + TMath::Sqrt((r-pt[2])*(r+pt[2]));
1091  c0 = (pt[0]*pt[0] + (pt[1]-r0)*(pt[1]+r0))/(dir[0]*dir[0]+dir[1]*dir[1]);
1092  delta = b0*b0-c0;
1093  if (delta>0) {
1094  y[nsol] = -b0-TMath::Sqrt(delta);
1095  if (y[nsol]>-tol) nsol++;
1096  y[nsol] = -b0+TMath::Sqrt(delta);
1097  if (y[nsol]>-tol) nsol++;
1098  }
1099  if (nsol) {
1100  // Sort solutions
1101  Int_t ind[4];
1102  TMath::Sort(nsol, y, ind,kFALSE);
1103  for (Int_t j=0; j<nsol; j++) x[j] = y[ind[j]];
1104  }
1105  } else {
1106  nsol = SolveQuartic(a,b,c,d,x);
1107  }
1108  if (!nsol) return TGeoShape::Big();
1109  // look for first positive solution
1110  Double_t phi, ndotd;
1111  Double_t r0[3], norm[3];
1112  Bool_t inner = (TMath::Abs(r-fRmin)<TGeoShape::Tolerance())?kTRUE:kFALSE;
1113  for (Int_t i=0; i<nsol; i++) {
1114  if (x[i]<-10) continue;
1115  phi = TMath::ATan2(pt[1]+x[i]*dir[1],pt[0]+x[i]*dir[0]);
1116  r0[0] = fR*TMath::Cos(phi);
1117  r0[1] = fR*TMath::Sin(phi);
1118  r0[2] = 0;
1119  for (Int_t ipt=0; ipt<3; ipt++) norm[ipt] = pt[ipt]+x[i]*dir[ipt] - r0[ipt];
1120  ndotd = norm[0]*dir[0]+norm[1]*dir[1]+norm[2]*dir[2];
1121  if (inner^in) {
1122  if (ndotd<0) continue;
1123  } else {
1124  if (ndotd>0) continue;
1125  }
1126  s = x[i];
1127  Double_t eps = TGeoShape::Big();
1128  Double_t delta = s*s*s*s + a*s*s*s + b*s*s + c*s + d;
1129  Double_t eps0 = -delta/(4.*s*s*s + 3.*a*s*s + 2.*b*s + c);
1130  while (TMath::Abs(eps)>TGeoShape::Tolerance()) {
1131  if (TMath::Abs(eps0)>100) break;
1132  s += eps0;
1133  if (TMath::Abs(s+eps0)<TGeoShape::Tolerance()) break;
1134  delta = s*s*s*s + a*s*s*s + b*s*s + c*s + d;
1135  eps = -delta/(4.*s*s*s + 3.*a*s*s + 2.*b*s + c);
1136  if (TMath::Abs(eps)>TMath::Abs(eps0)) break;
1137  eps0 = eps;
1138  }
1139  if (s<-TGeoShape::Tolerance()) continue;
1140  return TMath::Max(0.,s);
1141  }
1142  return TGeoShape::Big();
1143 }
1144 
1145 ////////////////////////////////////////////////////////////////////////////////
1146 /// Returns numbers of vertices, segments and polygons composing the shape mesh.
1147 
1148 void TGeoTorus::GetMeshNumbers(Int_t &nvert, Int_t &nsegs, Int_t &npols) const
1149 {
1150  Int_t n = gGeoManager->GetNsegments()+1;
1151  nvert = n*(n-1);
1152  Bool_t hasrmin = (GetRmin()>0)?kTRUE:kFALSE;
1153  Bool_t hasphi = (GetDphi()<360)?kTRUE:kFALSE;
1154  if (hasrmin) nvert *= 2;
1155  else if (hasphi) nvert += 2;
1156  nsegs = (2*n-1)*(n-1);
1157  npols = (n-1)*(n-1);
1158  if (hasrmin) {
1159  nsegs += (2*n-1)*(n-1);
1160  npols += (n-1)*(n-1);
1161  }
1162  if (hasphi) {
1163  nsegs += 2*(n-1);
1164  npols += 2*(n-1);
1165  }
1166 }
1167 
1168 ////////////////////////////////////////////////////////////////////////////////
1169 /// Fills a static 3D buffer and returns a reference.
1170 
1171 const TBuffer3D & TGeoTorus::GetBuffer3D(Int_t reqSections, Bool_t localFrame) const
1172 {
1173  static TBuffer3D buffer(TBuffer3DTypes::kGeneric);
1174 
1175  TGeoBBox::FillBuffer3D(buffer, reqSections, localFrame);
1176 
1177  if (reqSections & TBuffer3D::kRawSizes) {
1178  Int_t n = gGeoManager->GetNsegments()+1;
1179  Int_t nbPnts = n*(n-1);
1180  Bool_t hasrmin = (GetRmin()>0)?kTRUE:kFALSE;
1181  Bool_t hasphi = (GetDphi()<360)?kTRUE:kFALSE;
1182  if (hasrmin) nbPnts *= 2;
1183  else if (hasphi) nbPnts += 2;
1184 
1185  Int_t nbSegs = (2*n-1)*(n-1);
1186  Int_t nbPols = (n-1)*(n-1);
1187  if (hasrmin) {
1188  nbSegs += (2*n-1)*(n-1);
1189  nbPols += (n-1)*(n-1);
1190  }
1191  if (hasphi) {
1192  nbSegs += 2*(n-1);
1193  nbPols += 2*(n-1);
1194  }
1195 
1196  if (buffer.SetRawSizes(nbPnts, 3*nbPnts, nbSegs, 3*nbSegs, nbPols, 6*nbPols)) {
1197  buffer.SetSectionsValid(TBuffer3D::kRawSizes);
1198  }
1199  }
1200  // TODO: Push down to TGeoShape?? But would have to do raw sizes set first..
1201  // can rest of TGeoShape be deferred until after
1202  if ((reqSections & TBuffer3D::kRaw) && buffer.SectionsValid(TBuffer3D::kRawSizes)) {
1203  SetPoints(buffer.fPnts);
1204  if (!buffer.fLocalFrame) {
1205  TransformPoints(buffer.fPnts, buffer.NbPnts());
1206  }
1207 
1208  SetSegsAndPols(buffer);
1209  buffer.SetSectionsValid(TBuffer3D::kRaw);
1210  }
1211 
1212  return buffer;
1213 }
1214 
1215 ////////////////////////////////////////////////////////////////////////////////
1216 /// Check the inside status for each of the points in the array.
1217 /// Input: Array of point coordinates + vector size
1218 /// Output: Array of Booleans for the inside of each point
1219 
1220 void TGeoTorus::Contains_v(const Double_t *points, Bool_t *inside, Int_t vecsize) const
1221 {
1222  for (Int_t i=0; i<vecsize; i++) inside[i] = Contains(&points[3*i]);
1223 }
1224 
1225 ////////////////////////////////////////////////////////////////////////////////
1226 /// Compute the normal for an array o points so that norm.dot.dir is positive
1227 /// Input: Arrays of point coordinates and directions + vector size
1228 /// Output: Array of normal directions
1229 
1230 void TGeoTorus::ComputeNormal_v(const Double_t *points, const Double_t *dirs, Double_t *norms, Int_t vecsize)
1231 {
1232  for (Int_t i=0; i<vecsize; i++) ComputeNormal(&points[3*i], &dirs[3*i], &norms[3*i]);
1233 }
1234 
1235 ////////////////////////////////////////////////////////////////////////////////
1236 /// Compute distance from array of input points having directions specified by dirs. Store output in dists
1237 
1238 void TGeoTorus::DistFromInside_v(const Double_t *points, const Double_t *dirs, Double_t *dists, Int_t vecsize, Double_t* step) const
1239 {
1240  for (Int_t i=0; i<vecsize; i++) dists[i] = DistFromInside(&points[3*i], &dirs[3*i], 3, step[i]);
1241 }
1242 
1243 ////////////////////////////////////////////////////////////////////////////////
1244 /// Compute distance from array of input points having directions specified by dirs. Store output in dists
1245 
1246 void TGeoTorus::DistFromOutside_v(const Double_t *points, const Double_t *dirs, Double_t *dists, Int_t vecsize, Double_t* step) const
1247 {
1248  for (Int_t i=0; i<vecsize; i++) dists[i] = DistFromOutside(&points[3*i], &dirs[3*i], 3, step[i]);
1249 }
1250 
1251 ////////////////////////////////////////////////////////////////////////////////
1252 /// Compute safe distance from each of the points in the input array.
1253 /// Input: Array of point coordinates, array of statuses for these points, size of the arrays
1254 /// Output: Safety values
1255 
1256 void TGeoTorus::Safety_v(const Double_t *points, const Bool_t *inside, Double_t *safe, Int_t vecsize) const
1257 {
1258  for (Int_t i=0; i<vecsize; i++) safe[i] = Safety(&points[3*i], inside[i]);
1259 }