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