Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
TGeoParaboloid.cxx
Go to the documentation of this file.
1 // @(#)root/geom:$Id$
2 // Author: Mihaela Gheata 20/06/04
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 TGeoParaboloid
13 \ingroup Geometry_classes
14 
15 Paraboloid class.
16 
17 A paraboloid is the solid bounded by the following surfaces:
18  - 2 planes parallel with XY cutting the Z axis at Z=-dz and Z=+dz
19  - the surface of revolution of a parabola described by:
20  `z = a*(x*x + y*y) + b`
21 
22 The parameters a and b are automatically computed from:
23  - rlo - the radius of the circle of intersection between the
24  parabolic surface and the plane z = -dz
25  - rhi - the radius of the circle of intersection between the
26  parabolic surface and the plane z = +dz
27 
28  | -dz = a*rlo*rlo + b
29  | dz = a*rhi*rhi + b where: rlo != rhi, both >= 0
30 */
31 
32 #include "Riostream.h"
33 #include "TGeoManager.h"
34 #include "TGeoVolume.h"
35 #include "TVirtualGeoPainter.h"
36 #include "TGeoParaboloid.h"
37 #include "TVirtualPad.h"
38 #include "TBuffer3D.h"
39 #include "TBuffer3DTypes.h"
40 #include "TMath.h"
41 
42 ClassImp(TGeoParaboloid);
43 
44 ////////////////////////////////////////////////////////////////////////////////
45 /// Dummy constructor
46 
47 TGeoParaboloid::TGeoParaboloid()
48 {
49  fRlo = 0;
50  fRhi = 0;
51  fDz = 0;
52  fA = 0;
53  fB = 0;
54  SetShapeBit(TGeoShape::kGeoParaboloid);
55 }
56 
57 ////////////////////////////////////////////////////////////////////////////////
58 /// Default constructor specifying X and Y semiaxis length
59 
60 TGeoParaboloid::TGeoParaboloid(Double_t rlo, Double_t rhi, Double_t dz)
61  :TGeoBBox(0,0,0)
62 {
63  fRlo = 0;
64  fRhi = 0;
65  fDz = 0;
66  fA = 0;
67  fB = 0;
68  SetShapeBit(TGeoShape::kGeoParaboloid);
69  SetParaboloidDimensions(rlo, rhi, dz);
70  ComputeBBox();
71 }
72 
73 ////////////////////////////////////////////////////////////////////////////////
74 /// Default constructor specifying X and Y semiaxis length
75 
76 TGeoParaboloid::TGeoParaboloid(const char *name, Double_t rlo, Double_t rhi, Double_t dz)
77  :TGeoBBox(name, 0, 0, 0)
78 {
79  fRlo = 0;
80  fRhi = 0;
81  fDz = 0;
82  fA = 0;
83  fB = 0;
84  SetShapeBit(TGeoShape::kGeoParaboloid);
85  SetParaboloidDimensions(rlo, rhi, dz);
86  ComputeBBox();
87 }
88 
89 ////////////////////////////////////////////////////////////////////////////////
90 /// Default constructor specifying minimum and maximum radius
91 /// - param[0] = rlo
92 /// - param[1] = rhi
93 /// - param[2] = dz
94 
95 TGeoParaboloid::TGeoParaboloid(Double_t *param)
96 {
97  SetShapeBit(TGeoShape::kGeoParaboloid);
98  SetDimensions(param);
99  ComputeBBox();
100 }
101 
102 ////////////////////////////////////////////////////////////////////////////////
103 /// destructor
104 
105 TGeoParaboloid::~TGeoParaboloid()
106 {
107 }
108 
109 ////////////////////////////////////////////////////////////////////////////////
110 /// Computes capacity of the shape in [length^3]
111 
112 Double_t TGeoParaboloid::Capacity() const
113 {
114  Double_t capacity = TMath::Pi()*fDz*(fRlo*fRlo+fRhi*fRhi);
115  return capacity;
116 }
117 
118 ////////////////////////////////////////////////////////////////////////////////
119 /// compute bounding box of the tube
120 
121 void TGeoParaboloid::ComputeBBox()
122 {
123  fDX = TMath::Max(fRlo, fRhi);
124  fDY = fDX;
125  fDZ = fDz;
126 }
127 
128 ////////////////////////////////////////////////////////////////////////////////
129 /// Compute normal to closest surface from POINT.
130 
131 void TGeoParaboloid::ComputeNormal(const Double_t *point, const Double_t *dir, Double_t *norm)
132 {
133  norm[0] = norm[1] = 0.0;
134  if (TMath::Abs(point[2]) > fDz) {
135  norm[2] = TMath::Sign(1., dir[2]);
136  return;
137  }
138  Double_t safz = fDz-TMath::Abs(point[2]);
139  Double_t r = TMath::Sqrt(point[0]*point[0]+point[1]*point[1]);
140  Double_t safr = TMath::Abs(r-TMath::Sqrt((point[2]-fB)/fA));
141  if (safz<safr) {
142  norm[2] = TMath::Sign(1., dir[2]);
143  return;
144  }
145  Double_t talf = -2.*fA*r;
146  Double_t calf = 1./TMath::Sqrt(1.+talf*talf);
147  Double_t salf = talf * calf;
148  Double_t phi = TMath::ATan2(point[1], point[0]);
149 
150  norm[0] = salf*TMath::Cos(phi);
151  norm[1] = salf*TMath::Sin(phi);
152  norm[2] = calf;
153  Double_t ndotd = norm[0]*dir[0]+norm[1]*dir[1]+norm[2]*dir[2];
154  if (ndotd < 0) {
155  norm[0] = -norm[0];
156  norm[1] = -norm[1];
157  norm[2] = -norm[2];
158  }
159 }
160 
161 ////////////////////////////////////////////////////////////////////////////////
162 /// test if point is inside the elliptical tube
163 
164 Bool_t TGeoParaboloid::Contains(const Double_t *point) const
165 {
166  if (TMath::Abs(point[2])>fDz) return kFALSE;
167  Double_t aa = fA*(point[2]-fB);
168  if (aa < 0) return kFALSE;
169  Double_t rsq = point[0]*point[0]+point[1]*point[1];
170  if (aa < fA*fA*rsq) return kFALSE;
171  return kTRUE;
172 }
173 
174 ////////////////////////////////////////////////////////////////////////////////
175 /// compute closest distance from point px,py to each vertex
176 
177 Int_t TGeoParaboloid::DistancetoPrimitive(Int_t px, Int_t py)
178 {
179  Int_t n = gGeoManager->GetNsegments();
180  const Int_t numPoints=n*(n+1)+2;
181  return ShapeDistancetoPrimitive(numPoints, px, py);
182 }
183 
184 ////////////////////////////////////////////////////////////////////////////////
185 /// Compute distance from a point to the parabola given by:
186 /// `z = a*rsq + b; rsq = x*x+y*y`
187 
188 Double_t TGeoParaboloid::DistToParaboloid(const Double_t *point, const Double_t *dir, Bool_t in) const
189 {
190  Double_t rsq = point[0]*point[0]+point[1]*point[1];
191  Double_t a = fA * (dir[0]*dir[0] + dir[1]*dir[1]);
192  Double_t b = 2.*fA*(point[0]*dir[0]+point[1]*dir[1])-dir[2];
193  Double_t c = fA*rsq + fB - point[2];
194  Double_t dist = TGeoShape::Big();
195  if (TMath::Abs(a)<TGeoShape::Tolerance()) {
196  if (TMath::Abs(b)<TGeoShape::Tolerance()) return dist; // big
197  dist = -c/b;
198  if (dist < 0) return TGeoShape::Big();
199  return dist; // OK
200  }
201  Double_t ainv = 1./a;
202  Double_t sum = - b*ainv;
203  Double_t prod = c*ainv;
204  Double_t delta = sum*sum - 4.*prod;
205  if (delta<0) return dist; // big
206  delta = TMath::Sqrt(delta);
207  Double_t sone = TMath::Sign(1.,ainv);
208  Int_t i = -1;
209  while (i<2) {
210  dist = 0.5*(sum+i*sone*delta);
211  i += 2;
212  if (dist<0) continue;
213  if (dist<1.E-8) {
214  Double_t talf = -2.*fA*TMath::Sqrt(rsq);
215  Double_t phi = TMath::ATan2(point[1], point[0]);
216  Double_t ndotd = talf*(TMath::Cos(phi)*dir[0]+TMath::Sin(phi)*dir[1])+dir[2];
217  if (!in) ndotd *= -1;
218  if (ndotd<0) return dist;
219  } else return dist;
220  }
221  return TGeoShape::Big();
222 }
223 
224 ////////////////////////////////////////////////////////////////////////////////
225 /// compute distance from inside point to surface of the paraboloid
226 
227 Double_t TGeoParaboloid::DistFromInside(const Double_t *point, const Double_t *dir, Int_t iact, Double_t step, Double_t *safe) const
228 {
229  if (iact<3 && safe) {
230  // compute safe distance
231  *safe = Safety(point, kTRUE);
232  if (iact==0) return TGeoShape::Big();
233  if (iact==1 && step<*safe) return TGeoShape::Big();
234  }
235 
236  Double_t dz = TGeoShape::Big();
237  if (dir[2]<0) {
238  dz = -(point[2]+fDz)/dir[2];
239  } else if (dir[2]>0) {
240  dz = (fDz-point[2])/dir[2];
241  }
242  Double_t dpara = DistToParaboloid(point, dir, kTRUE);
243  return TMath::Min(dz, dpara);
244 }
245 
246 ////////////////////////////////////////////////////////////////////////////////
247 /// compute distance from outside point to surface of the paraboloid and safe distance
248 
249 Double_t TGeoParaboloid::DistFromOutside(const Double_t *point, const Double_t *dir, Int_t iact, Double_t step, Double_t *safe) const
250 {
251  Double_t snxt = TGeoShape::Big();
252  if (iact<3 && safe) {
253  // compute safe distance
254  *safe = Safety(point, kFALSE);
255  if (iact==0) return TGeoShape::Big();
256  if (iact==1 && step<*safe) return TGeoShape::Big();
257  }
258  Double_t xnew, ynew, znew;
259  if (point[2]<=-fDz) {
260  if (dir[2]<=0) return TGeoShape::Big();
261  snxt = -(fDz+point[2])/dir[2];
262  // find extrapolated X and Y
263  xnew = point[0]+snxt*dir[0];
264  ynew = point[1]+snxt*dir[1];
265  if ((xnew*xnew+ynew*ynew) <= fRlo*fRlo) return snxt;
266  } else if (point[2]>=fDz) {
267  if (dir[2]>=0) return TGeoShape::Big();
268  snxt = (fDz-point[2])/dir[2];
269  // find extrapolated X and Y
270  xnew = point[0]+snxt*dir[0];
271  ynew = point[1]+snxt*dir[1];
272  if ((xnew*xnew+ynew*ynew) <= fRhi*fRhi) return snxt;
273  }
274  snxt = DistToParaboloid(point, dir, kFALSE);
275  if (snxt > 1E20) return snxt;
276  znew = point[2]+snxt*dir[2];
277  if (TMath::Abs(znew) <= fDz) return snxt;
278  return TGeoShape::Big();
279 }
280 
281 ////////////////////////////////////////////////////////////////////////////////
282 /// Divide the paraboloid along one axis.
283 
284 TGeoVolume *TGeoParaboloid::Divide(TGeoVolume * /*voldiv*/, const char * /*divname*/, Int_t /*iaxis*/, Int_t /*ndiv*/,
285  Double_t /*start*/, Double_t /*step*/)
286 {
287  Error("Divide", "Paraboloid divisions not implemented");
288  return 0;
289 }
290 
291 ////////////////////////////////////////////////////////////////////////////////
292 /// Fill vector param[4] with the bounding cylinder parameters. The order
293 /// is the following : Rmin, Rmax, Phi1, Phi2
294 
295 void TGeoParaboloid::GetBoundingCylinder(Double_t *param) const
296 {
297  param[0] = 0.; // Rmin
298  param[1] = fDX; // Rmax
299  param[1] *= param[1];
300  param[2] = 0.; // Phi1
301  param[3] = 360.; // Phi2
302 }
303 
304 ////////////////////////////////////////////////////////////////////////////////
305 /// in case shape has some negative parameters, these has to be computed
306 /// in order to fit the mother
307 
308 TGeoShape *TGeoParaboloid::GetMakeRuntimeShape(TGeoShape *, TGeoMatrix *) const
309 {
310  return 0;
311 }
312 
313 ////////////////////////////////////////////////////////////////////////////////
314 /// print shape parameters
315 
316 void TGeoParaboloid::InspectShape() const
317 {
318  printf("*** Shape %s: TGeoParaboloid ***\n", GetName());
319  printf(" rlo = %11.5f\n", fRlo);
320  printf(" rhi = %11.5f\n", fRhi);
321  printf(" dz = %11.5f\n", fDz);
322  printf(" Bounding box:\n");
323  TGeoBBox::InspectShape();
324 }
325 
326 ////////////////////////////////////////////////////////////////////////////////
327 /// Creates a TBuffer3D describing *this* shape.
328 /// Coordinates are in local reference frame.
329 
330 TBuffer3D *TGeoParaboloid::MakeBuffer3D() const
331 {
332  Int_t n = gGeoManager->GetNsegments();
333  Int_t nbPnts = n*(n+1)+2;
334  Int_t nbSegs = n*(2*n+3);
335  Int_t nbPols = n*(n+2);
336 
337  TBuffer3D* buff = new TBuffer3D(TBuffer3DTypes::kGeneric,
338  nbPnts, 3*nbPnts, nbSegs, 3*nbSegs, nbPols, 2*n*5 + n*n*6);
339 
340  if (buff)
341  {
342  SetPoints(buff->fPnts);
343  SetSegsAndPols(*buff);
344  }
345 
346  return buff;
347 }
348 
349 ////////////////////////////////////////////////////////////////////////////////
350 /// Fill TBuffer3D structure for segments and polygons.
351 
352 void TGeoParaboloid::SetSegsAndPols(TBuffer3D &buff) const
353 {
354  Int_t indx, i, j;
355  Int_t n = gGeoManager->GetNsegments();
356 
357  Int_t c = GetBasicColor();
358 
359  Int_t nn1 = (n+1)*n+1;
360  indx = 0;
361  // Lower end-cap (n radial segments)
362  for (j=0; j<n; j++) {
363  buff.fSegs[indx++] = c+2;
364  buff.fSegs[indx++] = 0;
365  buff.fSegs[indx++] = j+1;
366  }
367  // Sectors (n)
368  for (i=0; i<n+1; i++) {
369  // lateral (circles) segments (n)
370  for (j=0; j<n; j++) {
371  buff.fSegs[indx++] = c;
372  buff.fSegs[indx++] = n*i+1+j;
373  buff.fSegs[indx++] = n*i+1+((j+1)%n);
374  }
375  if (i==n) break; // skip i=n for generators
376  // generator segments (n)
377  for (j=0; j<n; j++) {
378  buff.fSegs[indx++] = c;
379  buff.fSegs[indx++] = n*i+1+j;
380  buff.fSegs[indx++] = n*(i+1)+1+j;
381  }
382  }
383  // Upper end-cap
384  for (j=0; j<n; j++) {
385  buff.fSegs[indx++] = c+1;
386  buff.fSegs[indx++] = n*n+1+j;
387  buff.fSegs[indx++] = nn1;
388  }
389 
390  indx = 0;
391 
392  // lower end-cap (n polygons)
393  for (j=0; j<n; j++) {
394  buff.fPols[indx++] = c+2;
395  buff.fPols[indx++] = 3;
396  buff.fPols[indx++] = n+j;
397  buff.fPols[indx++] = (j+1)%n;
398  buff.fPols[indx++] = j;
399  }
400  // Sectors (n)
401  for (i=0; i<n; i++) {
402  // lateral faces (n)
403  for (j=0; j<n; j++) {
404  buff.fPols[indx++] = c;
405  buff.fPols[indx++] = 4;
406  buff.fPols[indx++] = (2*i+1)*n+j;
407  buff.fPols[indx++] = 2*(i+1)*n+j;
408  buff.fPols[indx++] = (2*i+3)*n+j;
409  buff.fPols[indx++] = 2*(i+1)*n+((j+1)%n);
410  }
411  }
412  // upper end-cap (n polygons)
413  for (j=0; j<n; j++) {
414  buff.fPols[indx++] = c+1;
415  buff.fPols[indx++] = 3;
416  buff.fPols[indx++] = 2*n*(n+1)+j;
417  buff.fPols[indx++] = 2*n*(n+1)+((j+1)%n);
418  buff.fPols[indx++] = (2*n+1)*n+j;
419  }
420 }
421 
422 ////////////////////////////////////////////////////////////////////////////////
423 /// Computes the closest distance from given point to this shape.
424 
425 Double_t TGeoParaboloid::Safety(const Double_t *point, Bool_t in) const
426 {
427  Double_t safz = fDz-TMath::Abs(point[2]);
428  if (!in) safz = -safz;
429  Double_t safr = TGeoShape::Big();
430  Double_t rsq = point[0]*point[0]+point[1]*point[1];
431  Double_t z0 = fA*rsq+fB;
432  Double_t r0sq = (point[2]-fB)/fA;
433  if (r0sq<0) {
434  if (in) return 0.;
435  return safz;
436  }
437  Double_t dr = TMath::Sqrt(rsq)-TMath::Sqrt(r0sq);
438  if (in) {
439  if (dr>-1.E-8) return 0.;
440  Double_t dz = TMath::Abs(point[2]-z0);
441  safr = -dr*dz/TMath::Sqrt(dr*dr+dz*dz);
442  } else {
443  if (dr<1.E-8) return safz;
444  Double_t talf = -2.*fA*TMath::Sqrt(r0sq);
445  Double_t salf = talf/TMath::Sqrt(1.+talf*talf);
446  safr = TMath::Abs(dr*salf);
447  }
448  if (in) return TMath::Min(safr,safz);
449  return TMath::Max(safr,safz);
450 }
451 
452 ////////////////////////////////////////////////////////////////////////////////
453 /// Set paraboloid dimensions.
454 
455 void TGeoParaboloid::SetParaboloidDimensions(Double_t rlo, Double_t rhi, Double_t dz)
456 {
457  if ((rlo<0) || (rhi<0) || (dz<=0) || TMath::Abs(rlo-rhi)<TGeoShape::Tolerance()) {
458  SetShapeBit(kGeoRunTimeShape);
459  Error("SetParaboloidDimensions", "Dimensions of %s invalid: check (rlo>=0) (rhi>=0) (rlo!=rhi) dz>0",GetName());
460  return;
461  }
462  fRlo = rlo;
463  fRhi = rhi;
464  fDz = dz;
465  Double_t dd = 1./(fRhi*fRhi - fRlo*fRlo);
466  fA = 2.*fDz*dd;
467  fB = - fDz * (fRlo*fRlo + fRhi*fRhi)*dd;
468 }
469 
470 ////////////////////////////////////////////////////////////////////////////////
471 /// Set paraboloid dimensions starting from an array.
472 
473 void TGeoParaboloid::SetDimensions(Double_t *param)
474 {
475  Double_t rlo = param[0];
476  Double_t rhi = param[1];
477  Double_t dz = param[2];
478  SetParaboloidDimensions(rlo, rhi, dz);
479 }
480 
481 ////////////////////////////////////////////////////////////////////////////////
482 /// Create paraboloid mesh points.
483 /// ~~~ {.cpp}
484 /// Npoints = n*(n+1) + 2
485 /// ifirst = 0
486 /// ipoint(i,j) = 1+i*n+j; i=[0,n] j=[0,n-1]
487 /// ilast = 1+n*(n+1)
488 /// Nsegments = n*(2*n+3)
489 /// lower: (0, j+1); j=[0,n-1]
490 /// circle(i): (n*i+1+j, n*i+1+(j+1)%n); i=[0,n] j=[0,n-1]
491 /// generator(i): (n*i+1+j, n*(i+1)+1+j); i,j=[0,n-1]
492 /// upper: (n*n+1+j, (n+1)*n+1) j=[0,n-1]
493 /// Npolygons = n*(n+2)
494 /// lower: (n+j, (j+1)%n, j) j=[0,n-1]
495 /// lateral(i): ((2*i+1)*n+j, 2*(i+1)*n+j, (2*i+3)*n+j, 2*(i+1)*n+(j+1)%n)
496 /// i,j = [0,n-1]
497 /// upper: ((2n+1)*n+j, 2*n*(n+1)+(j+1)%n, 2*n*(n+1)+j) j=[0,n-1]
498 /// ~~~
499 
500 void TGeoParaboloid::SetPoints(Double_t *points) const
501 {
502  if (!points) return;
503  Double_t ttmin, ttmax;
504  ttmin = TMath::ATan2(-fDz, fRlo);
505  ttmax = TMath::ATan2(fDz, fRhi);
506  Int_t n = gGeoManager->GetNsegments();
507  Double_t dtt = (ttmax-ttmin)/n;
508  Double_t dphi = 360./n;
509  Double_t tt;
510  Double_t r, z, delta;
511  Double_t phi, sph, cph;
512  Int_t indx = 0;
513  // center of the lower endcap:
514  points[indx++] = 0; // x
515  points[indx++] = 0; // y
516  points[indx++] = -fDz;
517  for (Int_t i=0; i<n+1; i++) { // nz planes = n+1
518  if (i==0) {
519  r = fRlo;
520  z = -fDz;
521  } else if (i==n) {
522  r = fRhi;
523  z = fDz;
524  } else {
525  tt = TMath::Tan(ttmin + i*dtt);
526  delta = tt*tt - 4*fA*fB; // should be always positive (a*b<0)
527  r = 0.5*(tt+TMath::Sqrt(delta))/fA;
528  z = r*tt;
529  }
530  for (Int_t j=0; j<n; j++) {
531  phi = j*dphi*TMath::DegToRad();
532  sph=TMath::Sin(phi);
533  cph=TMath::Cos(phi);
534  points[indx++] = r*cph;
535  points[indx++] = r*sph;
536  points[indx++] = z;
537  }
538  }
539  // center of the upper endcap
540  points[indx++] = 0; // x
541  points[indx++] = 0; // y
542  points[indx++] = fDz;
543 }
544 
545 ////////////////////////////////////////////////////////////////////////////////
546 /// Returns numbers of vertices, segments and polygons composing the shape mesh.
547 
548 void TGeoParaboloid::GetMeshNumbers(Int_t &nvert, Int_t &nsegs, Int_t &npols) const
549 {
550  Int_t n = gGeoManager->GetNsegments();
551  nvert = n*(n+1)+2;
552  nsegs = n*(2*n+3);
553  npols = n*(n+2);
554 }
555 
556 ////////////////////////////////////////////////////////////////////////////////
557 /// Returns number of vertices on the paraboloid mesh.
558 
559 Int_t TGeoParaboloid::GetNmeshVertices() const
560 {
561  Int_t n = gGeoManager->GetNsegments();
562  return (n*(n+1)+2);
563 }
564 
565 ////////////////////////////////////////////////////////////////////////////////
566 /// Save a primitive as a C++ statement(s) on output stream "out".
567 
568 void TGeoParaboloid::SavePrimitive(std::ostream &out, Option_t * /*option*/ /*= ""*/)
569 {
570  if (TObject::TestBit(kGeoSavePrimitive)) return;
571  out << " // Shape: " << GetName() << " type: " << ClassName() << std::endl;
572  out << " rlo = " << fRlo << ";" << std::endl;
573  out << " rhi = " << fRhi << ";" << std::endl;
574  out << " dz = " << fDZ << ";" << std::endl;
575  out << " TGeoShape *" << GetPointerName() << " = new TGeoParaboloid(\"" << GetName() << "\", rlo,rhi,dz);" << std::endl;
576  TObject::SetBit(TGeoShape::kGeoSavePrimitive);
577 }
578 
579 ////////////////////////////////////////////////////////////////////////////////
580 /// Create paraboloid mesh points.
581 
582 void TGeoParaboloid::SetPoints(Float_t *points) const
583 {
584  if (!points) return;
585  Double_t ttmin, ttmax;
586  ttmin = TMath::ATan2(-fDz, fRlo);
587  ttmax = TMath::ATan2(fDz, fRhi);
588  Int_t n = gGeoManager->GetNsegments();
589  Double_t dtt = (ttmax-ttmin)/n;
590  Double_t dphi = 360./n;
591  Double_t tt;
592  Double_t r, z, delta;
593  Double_t phi, sph, cph;
594  Int_t indx = 0;
595  // center of the lower endcap:
596  points[indx++] = 0; // x
597  points[indx++] = 0; // y
598  points[indx++] = -fDz;
599  for (Int_t i=0; i<n+1; i++) { // nz planes = n+1
600  if (i==0) {
601  r = fRlo;
602  z = -fDz;
603  } else if (i==n) {
604  r = fRhi;
605  z = fDz;
606  } else {
607  tt = TMath::Tan(ttmin + i*dtt);
608  delta = tt*tt - 4*fA*fB; // should be always positive (a*b<0)
609  r = 0.5*(tt+TMath::Sqrt(delta))/fA;
610  z = r*tt;
611  }
612  for (Int_t j=0; j<n; j++) {
613  phi = j*dphi*TMath::DegToRad();
614  sph=TMath::Sin(phi);
615  cph=TMath::Cos(phi);
616  points[indx++] = r*cph;
617  points[indx++] = r*sph;
618  points[indx++] = z;
619  }
620  }
621  // center of the upper endcap
622  points[indx++] = 0; // x
623  points[indx++] = 0; // y
624  points[indx++] = fDz;
625 }
626 
627 ////////////////////////////////////////////////////////////////////////////////
628 void TGeoParaboloid::Sizeof3D() const
629 {
630 }
631 
632 ////////////////////////////////////////////////////////////////////////////////
633 /// Fills a static 3D buffer and returns a reference.
634 
635 const TBuffer3D & TGeoParaboloid::GetBuffer3D(Int_t reqSections, Bool_t localFrame) const
636 {
637  static TBuffer3D buffer(TBuffer3DTypes::kGeneric);
638  TGeoBBox::FillBuffer3D(buffer, reqSections, localFrame);
639 
640  if (reqSections & TBuffer3D::kRawSizes) {
641  Int_t n = gGeoManager->GetNsegments();
642  Int_t nbPnts = n*(n+1)+2;
643  Int_t nbSegs = n*(2*n+3);
644  Int_t nbPols = n*(n+2);
645  if (buffer.SetRawSizes(nbPnts, 3*nbPnts, nbSegs, 3*nbSegs, nbPols, 2*n*5 + n*n*6)) {
646  buffer.SetSectionsValid(TBuffer3D::kRawSizes);
647  }
648  }
649  if ((reqSections & TBuffer3D::kRaw) && buffer.SectionsValid(TBuffer3D::kRawSizes)) {
650  SetPoints(buffer.fPnts);
651  if (!buffer.fLocalFrame) {
652  TransformPoints(buffer.fPnts, buffer.NbPnts());
653  }
654  SetSegsAndPols(buffer);
655  buffer.SetSectionsValid(TBuffer3D::kRaw);
656  }
657 
658  return buffer;
659 }
660 
661 ////////////////////////////////////////////////////////////////////////////////
662 /// Check the inside status for each of the points in the array.
663 /// Input: Array of point coordinates + vector size
664 /// Output: Array of Booleans for the inside of each point
665 
666 void TGeoParaboloid::Contains_v(const Double_t *points, Bool_t *inside, Int_t vecsize) const
667 {
668  for (Int_t i=0; i<vecsize; i++) inside[i] = Contains(&points[3*i]);
669 }
670 
671 ////////////////////////////////////////////////////////////////////////////////
672 /// Compute the normal for an array o points so that norm.dot.dir is positive
673 /// Input: Arrays of point coordinates and directions + vector size
674 /// Output: Array of normal directions
675 
676 void TGeoParaboloid::ComputeNormal_v(const Double_t *points, const Double_t *dirs, Double_t *norms, Int_t vecsize)
677 {
678  for (Int_t i=0; i<vecsize; i++) ComputeNormal(&points[3*i], &dirs[3*i], &norms[3*i]);
679 }
680 
681 ////////////////////////////////////////////////////////////////////////////////
682 /// Compute distance from array of input points having directions specified by dirs. Store output in dists
683 
684 void TGeoParaboloid::DistFromInside_v(const Double_t *points, const Double_t *dirs, Double_t *dists, Int_t vecsize, Double_t* step) const
685 {
686  for (Int_t i=0; i<vecsize; i++) dists[i] = DistFromInside(&points[3*i], &dirs[3*i], 3, step[i]);
687 }
688 
689 ////////////////////////////////////////////////////////////////////////////////
690 /// Compute distance from array of input points having directions specified by dirs. Store output in dists
691 
692 void TGeoParaboloid::DistFromOutside_v(const Double_t *points, const Double_t *dirs, Double_t *dists, Int_t vecsize, Double_t* step) const
693 {
694  for (Int_t i=0; i<vecsize; i++) dists[i] = DistFromOutside(&points[3*i], &dirs[3*i], 3, step[i]);
695 }
696 
697 ////////////////////////////////////////////////////////////////////////////////
698 /// Compute safe distance from each of the points in the input array.
699 /// Input: Array of point coordinates, array of statuses for these points, size of the arrays
700 /// Output: Safety values
701 
702 void TGeoParaboloid::Safety_v(const Double_t *points, const Bool_t *inside, Double_t *safe, Int_t vecsize) const
703 {
704  for (Int_t i=0; i<vecsize; i++) safe[i] = Safety(&points[3*i], inside[i]);
705 }