Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
TGeoCone.cxx
Go to the documentation of this file.
1 // @(#)root/geom:$Id$
2 // Author: Andrei Gheata 31/01/02
3 // TGeoCone::Contains() and DistFromInside() implemented by Mihaela Gheata
4 
5 /*************************************************************************
6  * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers. *
7  * All rights reserved. *
8  * *
9  * For the licensing terms see $ROOTSYS/LICENSE. *
10  * For the list of contributors see $ROOTSYS/README/CREDITS. *
11  *************************************************************************/
12 
13 /** \class TGeoCone
14 \ingroup Geometry_classes
15 
16 Conical tube class. It has 5 parameters :
17  - dz - half length in z
18  - Rmin1, Rmax1 - inside and outside radii at -dz
19  - Rmin2, Rmax2 - inside and outside radii at +dz
20 
21 Begin_Macro(source)
22 {
23  TCanvas *c = new TCanvas("c", "c",0,0,600,600);
24  new TGeoManager("cone", "poza4");
25  TGeoMaterial *mat = new TGeoMaterial("Al", 26.98,13,2.7);
26  TGeoMedium *med = new TGeoMedium("MED",1,mat);
27  TGeoVolume *top = gGeoManager->MakeBox("TOP",med,100,100,100);
28  gGeoManager->SetTopVolume(top);
29  TGeoVolume *vol = gGeoManager->MakeCone("CONE",med, 40,10,20,35,45);
30  vol->SetLineWidth(2);
31  top->AddNode(vol,1);
32  gGeoManager->CloseGeometry();
33  gGeoManager->SetNsegments(30);
34  top->Draw();
35  TView *view = gPad->GetView();
36  view->ShowAxis();
37 }
38 End_Macro
39 */
40 
41 
42 /** \class TGeoConeSeg
43 \ingroup Geometry_classes
44 
45 A phi segment of a conical tube. Has 7 parameters :
46  - the same 5 as a cone;
47  - first phi limit (in degrees)
48  - second phi limit
49 
50 Begin_Macro(source)
51 {
52  TCanvas *c = new TCanvas("c", "c",0,0,600,600);
53  new TGeoManager("coneseg", "poza5");
54  TGeoMaterial *mat = new TGeoMaterial("Al", 26.98,13,2.7);
55  TGeoMedium *med = new TGeoMedium("MED",1,mat);
56  TGeoVolume *top = gGeoManager->MakeBox("TOP",med,100,100,100);
57  gGeoManager->SetTopVolume(top);
58  TGeoVolume *vol = gGeoManager->MakeCons("CONESEG",med, 40,30,40,10,20,-30,250);
59  top->AddNode(vol,1);
60  gGeoManager->CloseGeometry();
61  gGeoManager->SetNsegments(30);
62  top->Draw();
63  TView *view = gPad->GetView();
64  view->ShowAxis();
65 }
66 End_Macro
67 */
68 
69 #include "Riostream.h"
70 
71 #include "TGeoManager.h"
72 #include "TGeoVolume.h"
73 #include "TVirtualGeoPainter.h"
74 #include "TGeoCone.h"
75 #include "TVirtualPad.h"
76 #include "TBuffer3D.h"
77 #include "TBuffer3DTypes.h"
78 #include "TMath.h"
79 
80 ClassImp(TGeoCone);
81 
82 ////////////////////////////////////////////////////////////////////////////////
83 /// Default constructor
84 
85 TGeoCone::TGeoCone()
86 {
87  SetShapeBit(TGeoShape::kGeoCone);
88  fDz = 0.0;
89  fRmin1 = 0.0;
90  fRmax1 = 0.0;
91  fRmin2 = 0.0;
92  fRmax2 = 0.0;
93 }
94 
95 ////////////////////////////////////////////////////////////////////////////////
96 /// Default constructor specifying minimum and maximum radius
97 
98 TGeoCone::TGeoCone(Double_t dz, Double_t rmin1, Double_t rmax1,
99  Double_t rmin2, Double_t rmax2)
100  :TGeoBBox(0, 0, 0)
101 {
102  SetShapeBit(TGeoShape::kGeoCone);
103  SetConeDimensions(dz, rmin1, rmax1, rmin2, rmax2);
104  if ((dz<0) || (rmin1<0) || (rmax1<0) || (rmin2<0) || (rmax2<0)) {
105  SetShapeBit(kGeoRunTimeShape);
106  }
107  else ComputeBBox();
108 }
109 
110 ////////////////////////////////////////////////////////////////////////////////
111 /// Default constructor specifying minimum and maximum radius
112 
113 TGeoCone::TGeoCone(const char *name, Double_t dz, Double_t rmin1, Double_t rmax1,
114  Double_t rmin2, Double_t rmax2)
115  :TGeoBBox(name, 0, 0, 0)
116 {
117  SetShapeBit(TGeoShape::kGeoCone);
118  SetConeDimensions(dz, rmin1, rmax1, rmin2, rmax2);
119  if ((dz<0) || (rmin1<0) || (rmax1<0) || (rmin2<0) || (rmax2<0)) {
120  SetShapeBit(kGeoRunTimeShape);
121  }
122  else ComputeBBox();
123 }
124 
125 ////////////////////////////////////////////////////////////////////////////////
126 /// Default constructor specifying minimum and maximum radius
127 /// - param[0] = dz
128 /// - param[1] = Rmin1
129 /// - param[2] = Rmax1
130 /// - param[3] = Rmin2
131 /// - param[4] = Rmax2
132 
133 TGeoCone::TGeoCone(Double_t *param)
134  :TGeoBBox(0, 0, 0)
135 {
136  SetShapeBit(TGeoShape::kGeoCone);
137  SetDimensions(param);
138  if ((fDz<0) || (fRmin1<0) || (fRmax1<0) || (fRmin2<0) || (fRmax2<0))
139  SetShapeBit(kGeoRunTimeShape);
140  else ComputeBBox();
141 }
142 
143 ////////////////////////////////////////////////////////////////////////////////
144 /// Computes capacity of the shape in [length^3]
145 
146 Double_t TGeoCone::Capacity() const
147 {
148  return TGeoCone::Capacity(fDz, fRmin1, fRmax1, fRmin2, fRmax2);
149 }
150 
151 ////////////////////////////////////////////////////////////////////////////////
152 /// Computes capacity of the shape in [length^3]
153 
154 Double_t TGeoCone::Capacity(Double_t dz, Double_t rmin1, Double_t rmax1, Double_t rmin2, Double_t rmax2)
155 {
156  Double_t capacity = (2.*dz*TMath::Pi()/3.)*(rmax1*rmax1+rmax2*rmax2+rmax1*rmax2-
157  rmin1*rmin1-rmin2*rmin2-rmin1*rmin2);
158  return capacity;
159 }
160 
161 ////////////////////////////////////////////////////////////////////////////////
162 /// destructor
163 
164 TGeoCone::~TGeoCone()
165 {
166 }
167 
168 ////////////////////////////////////////////////////////////////////////////////
169 /// compute bounding box of the sphere
170 
171 void TGeoCone::ComputeBBox()
172 {
173  TGeoBBox *box = (TGeoBBox*)this;
174  box->SetBoxDimensions(TMath::Max(fRmax1, fRmax2), TMath::Max(fRmax1, fRmax2), fDz);
175  memset(fOrigin, 0, 3*sizeof(Double_t));
176 }
177 
178 ////////////////////////////////////////////////////////////////////////////////
179 /// Compute normal to closest surface from POINT.
180 
181 void TGeoCone::ComputeNormal(const Double_t *point, const Double_t *dir, Double_t *norm)
182 {
183  Double_t safr,safe,phi;
184  memset(norm,0,3*sizeof(Double_t));
185  phi = TMath::ATan2(point[1],point[0]);
186  Double_t cphi = TMath::Cos(phi);
187  Double_t sphi = TMath::Sin(phi);
188  Double_t ro1 = 0.5*(fRmin1+fRmin2);
189  Double_t tg1 = 0.5*(fRmin2-fRmin1)/fDz;
190  Double_t cr1 = 1./TMath::Sqrt(1.+tg1*tg1);
191  Double_t ro2 = 0.5*(fRmax1+fRmax2);
192  Double_t tg2 = 0.5*(fRmax2-fRmax1)/fDz;
193  Double_t cr2 = 1./TMath::Sqrt(1.+tg2*tg2);
194 
195  Double_t r=TMath::Sqrt(point[0]*point[0]+point[1]*point[1]);
196  Double_t rin = tg1*point[2]+ro1;
197  Double_t rout = tg2*point[2]+ro2;
198  safe = TMath::Abs(fDz-TMath::Abs(point[2]));
199  norm[2] = 1;
200 
201  safr = (ro1>0)?(TMath::Abs((r-rin)*cr1)):TGeoShape::Big();
202  if (safr<safe) {
203  safe = safr;
204  norm[0] = cr1*cphi;
205  norm[1] = cr1*sphi;
206  norm[2] = -tg1*cr1;
207  }
208  safr = TMath::Abs((rout-r)*cr2);
209  if (safr<safe) {
210  norm[0] = cr2*cphi;
211  norm[1] = cr2*sphi;
212  norm[2] = -tg2*cr2;
213  }
214  if (norm[0]*dir[0]+norm[1]*dir[1]+norm[2]*dir[2]<0) {
215  norm[0] = -norm[0];
216  norm[1] = -norm[1];
217  norm[2] = -norm[2];
218  }
219 }
220 
221 ////////////////////////////////////////////////////////////////////////////////
222 /// Compute normal to closest surface from POINT.
223 
224 void TGeoCone::ComputeNormalS(const Double_t *point, const Double_t *dir, Double_t *norm,
225  Double_t dz, Double_t rmin1, Double_t rmax1, Double_t rmin2, Double_t rmax2)
226 {
227  Double_t safe,phi;
228  memset(norm,0,3*sizeof(Double_t));
229  phi = TMath::ATan2(point[1],point[0]);
230  Double_t cphi = TMath::Cos(phi);
231  Double_t sphi = TMath::Sin(phi);
232  Double_t ro1 = 0.5*(rmin1+rmin2);
233  Double_t tg1 = 0.5*(rmin2-rmin1)/dz;
234  Double_t cr1 = 1./TMath::Sqrt(1.+tg1*tg1);
235  Double_t ro2 = 0.5*(rmax1+rmax2);
236  Double_t tg2 = 0.5*(rmax2-rmax1)/dz;
237  Double_t cr2 = 1./TMath::Sqrt(1.+tg2*tg2);
238 
239  Double_t r=TMath::Sqrt(point[0]*point[0]+point[1]*point[1]);
240  Double_t rin = tg1*point[2]+ro1;
241  Double_t rout = tg2*point[2]+ro2;
242  safe = (ro1>0)?(TMath::Abs((r-rin)*cr1)):TGeoShape::Big();
243  norm[0] = cr1*cphi;
244  norm[1] = cr1*sphi;
245  norm[2] = -tg1*cr1;
246  if (TMath::Abs((rout-r)*cr2)<safe) {
247  norm[0] = cr2*cphi;
248  norm[1] = cr2*sphi;
249  norm[2] = -tg2*cr2;
250  }
251  if (norm[0]*dir[0]+norm[1]*dir[1]+norm[2]*dir[2]<0) {
252  norm[0] = -norm[0];
253  norm[1] = -norm[1];
254  norm[2] = -norm[2];
255  }
256 }
257 
258 ////////////////////////////////////////////////////////////////////////////////
259 /// test if point is inside this cone
260 
261 Bool_t TGeoCone::Contains(const Double_t *point) const
262 {
263  if (TMath::Abs(point[2]) > fDz) return kFALSE;
264  Double_t r2 = point[0]*point[0]+point[1]*point[1];
265  Double_t rl = 0.5*(fRmin2*(point[2]+fDz)+fRmin1*(fDz-point[2]))/fDz;
266  Double_t rh = 0.5*(fRmax2*(point[2]+fDz)+fRmax1*(fDz-point[2]))/fDz;
267  if ((r2<rl*rl) || (r2>rh*rh)) return kFALSE;
268  return kTRUE;
269 }
270 
271 ////////////////////////////////////////////////////////////////////////////////
272 /// Compute distance from inside point to surface of the cone (static)
273 /// Boundary safe algorithm.
274 
275 Double_t TGeoCone::DistFromInsideS(const Double_t *point, const Double_t *dir, Double_t dz,
276  Double_t rmin1, Double_t rmax1, Double_t rmin2, Double_t rmax2)
277 {
278  if (dz<=0) return TGeoShape::Big();
279  // compute distance to surface
280  // Do Z
281  Double_t sz = TGeoShape::Big();
282  if (dir[2]) {
283  sz = (TMath::Sign(dz, dir[2])-point[2])/dir[2];
284  if (sz<=0) return 0.0;
285  }
286  Double_t rsq=point[0]*point[0]+point[1]*point[1];
287  Double_t zinv = 1./dz;
288  Double_t rin = 0.5*(rmin1+rmin2+(rmin2-rmin1)*point[2]*zinv);
289  // Do Rmin
290  Double_t sr = TGeoShape::Big();
291  Double_t b,delta,zi;
292  if (rin>0) {
293  // Protection in case point is actually outside the cone
294  if (rsq < rin*(rin+TGeoShape::Tolerance())) {
295  Double_t ddotn = point[0]*dir[0]+point[1]*dir[1]+0.5*(rmin1-rmin2)*dir[2]*zinv*TMath::Sqrt(rsq);
296  if (ddotn<=0) return 0.0;
297  } else {
298  TGeoCone::DistToCone(point, dir, dz, rmin1, rmin2, b, delta);
299  if (delta>0) {
300  sr = -b-delta;
301  if (sr>0) {
302  zi = point[2]+sr*dir[2];
303  if (TMath::Abs(zi)<=dz) return TMath::Min(sz,sr);
304  }
305  sr = -b+delta;
306  if (sr>0) {
307  zi = point[2]+sr*dir[2];
308  if (TMath::Abs(zi)<=dz) return TMath::Min(sz,sr);
309  }
310  }
311  }
312  }
313  // Do Rmax
314  Double_t rout = 0.5*(rmax1+rmax2+(rmax2-rmax1)*point[2]*zinv);
315  if (rsq > rout*(rout-TGeoShape::Tolerance())) {
316  Double_t ddotn = point[0]*dir[0]+point[1]*dir[1]+0.5*(rmax1-rmax2)*dir[2]*zinv*TMath::Sqrt(rsq);
317  if (ddotn>=0) return 0.0;
318  TGeoCone::DistToCone(point, dir, dz, rmax1, rmax2, b, delta);
319  if (delta<0) return 0.0;
320  sr = -b+delta;
321  if (sr<0) return sz;
322  if (TMath::Abs(-b-delta)>sr) return sz;
323  zi = point[2]+sr*dir[2];
324  if (TMath::Abs(zi)<=dz) return TMath::Min(sz,sr);
325  return sz;
326  }
327  TGeoCone::DistToCone(point, dir, dz, rmax1, rmax2, b, delta);
328  if (delta>0) {
329  sr = -b-delta;
330  if (sr>0) {
331  zi = point[2]+sr*dir[2];
332  if (TMath::Abs(zi)<=dz) return TMath::Min(sz,sr);
333  }
334  sr = -b+delta;
335  if (sr>TGeoShape::Tolerance()) {
336  zi = point[2]+sr*dir[2];
337  if (TMath::Abs(zi)<=dz) return TMath::Min(sz,sr);
338  }
339  }
340  return sz;
341 }
342 
343 ////////////////////////////////////////////////////////////////////////////////
344 /// Compute distance from inside point to surface of the cone
345 /// Boundary safe algorithm.
346 
347 Double_t TGeoCone::DistFromInside(const Double_t *point, const Double_t *dir, Int_t iact, Double_t step, Double_t *safe) const
348 {
349  if (iact<3 && safe) {
350  *safe = Safety(point, kTRUE);
351  if (iact==0) return TGeoShape::Big();
352  if ((iact==1) && (*safe>step)) return TGeoShape::Big();
353  }
354  // compute distance to surface
355  return TGeoCone::DistFromInsideS(point, dir, fDz, fRmin1, fRmax1, fRmin2, fRmax2);
356 }
357 
358 ////////////////////////////////////////////////////////////////////////////////
359 /// Compute distance from outside point to surface of the tube
360 /// Boundary safe algorithm.
361 
362 Double_t TGeoCone::DistFromOutsideS(const Double_t *point, const Double_t *dir, Double_t dz,
363  Double_t rmin1, Double_t rmax1, Double_t rmin2, Double_t rmax2)
364 {
365  // compute distance to Z planes
366  if (dz<=0) return TGeoShape::Big();
367  Double_t snxt;
368  Double_t xp, yp, zp;
369  Bool_t inz = kTRUE;
370 
371  if (point[2]<=-dz) {
372  if (dir[2]<=0) return TGeoShape::Big();
373  snxt = (-dz-point[2])/dir[2];
374  xp = point[0]+snxt*dir[0];
375  yp = point[1]+snxt*dir[1];
376  Double_t r2 = xp*xp+yp*yp;
377  if ((r2>=rmin1*rmin1) && (r2<=rmax1*rmax1)) return snxt;
378  inz = kFALSE;
379  } else {
380  if (point[2]>=dz) {
381  if (dir[2]>=0) return TGeoShape::Big();
382  snxt = (dz-point[2])/dir[2];
383  xp = point[0]+snxt*dir[0];
384  yp = point[1]+snxt*dir[1];
385  Double_t r2 = xp*xp+yp*yp;
386  if ((r2>=rmin2*rmin2) && (r2<=rmax2*rmax2)) return snxt;
387  inz = kFALSE;
388  }
389  }
390 
391  Double_t rsq = point[0]*point[0]+point[1]*point[1];
392  Double_t dzinv = 1./dz;
393  Double_t ro1=0.5*(rmin1+rmin2);
394  Bool_t hasrmin = (ro1>0)?kTRUE:kFALSE;
395  Double_t tg1 = 0.;
396  Double_t rin = 0.;
397  Bool_t inrmin = kTRUE; // r>=rmin
398  if (hasrmin) {
399  tg1=0.5*(rmin2-rmin1)*dzinv;
400  rin=ro1+tg1*point[2];
401  if (rin>0 && rsq<rin*(rin-TGeoShape::Tolerance())) inrmin=kFALSE;
402  }
403  Double_t ro2=0.5*(rmax1+rmax2);
404  Double_t tg2=0.5*(rmax2-rmax1)*dzinv;
405  Double_t rout=tg2*point[2]+ro2;
406  Bool_t inrmax = kFALSE;
407  if (rout>0 && rsq<rout*(rout+TGeoShape::Tolerance())) inrmax=kTRUE;
408  Bool_t in = inz & inrmin & inrmax;
409  Double_t b,delta;
410  // If inside cone, we are most likely on a boundary within machine precision.
411  if (in) {
412  Double_t r=TMath::Sqrt(rsq);
413  Double_t safz = dz-TMath::Abs(point[2]); // positive
414  Double_t safrmin = (hasrmin)?(r-rin):TGeoShape::Big();
415  Double_t safrmax = rout - r;
416  if (safz<=safrmin && safz<=safrmax) {
417  // on Z boundary
418  if (point[2]*dir[2]<0) return 0.0;
419  return TGeoShape::Big();
420  }
421  if (safrmax<safrmin) {
422  // on rmax boundary
423  Double_t ddotn = point[0]*dir[0]+point[1]*dir[1]-tg2*dir[2]*r;
424  if (ddotn<=0) return 0.0;
425  return TGeoShape::Big();
426  }
427  // on rmin boundary
428  Double_t ddotn = point[0]*dir[0]+point[1]*dir[1]-tg1*dir[2]*r;
429  if (ddotn>=0) return 0.0;
430  // we can cross (+) solution of rmin
431  TGeoCone::DistToCone(point, dir, dz, rmin1, rmin2, b, delta);
432 
433  if (delta<0) return 0.0;
434  snxt = -b+delta;
435  if (snxt<0) return TGeoShape::Big();
436  if (TMath::Abs(-b-delta)>snxt) return TGeoShape::Big();
437  zp = point[2]+snxt*dir[2];
438  if (TMath::Abs(zp)<=dz) return snxt;
439  return TGeoShape::Big();
440  }
441 
442  // compute distance to inner cone
443  snxt = TGeoShape::Big();
444  if (!inrmin) {
445  // ray can cross inner cone (but not only!)
446  TGeoCone::DistToCone(point, dir, dz, rmin1, rmin2, b, delta);
447  if (delta<0) return TGeoShape::Big();
448  snxt = -b+delta;
449  if (snxt>0) {
450  zp = point[2]+snxt*dir[2];
451  if (TMath::Abs(zp)<=dz) return snxt;
452  }
453  snxt = -b-delta;
454  if (snxt>0) {
455  zp = point[2]+snxt*dir[2];
456  if (TMath::Abs(zp)<=dz) return snxt;
457  }
458  snxt = TGeoShape::Big();
459  } else {
460  if (hasrmin) {
461  TGeoCone::DistToCone(point, dir, dz, rmin1, rmin2, b, delta);
462  if (delta>0) {
463  Double_t din = -b+delta;
464  if (din>0) {
465  zp = point[2]+din*dir[2];
466  if (TMath::Abs(zp)<=dz) snxt = din;
467  }
468  }
469  }
470  }
471 
472  if (inrmax) return snxt;
473  // We can cross outer cone, both solutions possible
474  // compute distance to outer cone
475  TGeoCone::DistToCone(point, dir, dz, rmax1, rmax2, b, delta);
476  if (delta<0) return snxt;
477  Double_t dout = -b-delta;
478  if (dout>0 && dout<snxt) {
479  zp = point[2]+dout*dir[2];
480  if (TMath::Abs(zp)<=dz) return dout;
481  }
482  dout = -b+delta;
483  if (dout<=0 || dout>snxt) return snxt;
484  zp = point[2]+dout*dir[2];
485  if (TMath::Abs(zp)<=dz) return dout;
486  return snxt;
487 }
488 
489 ////////////////////////////////////////////////////////////////////////////////
490 /// compute distance from outside point to surface of the tube
491 
492 Double_t TGeoCone::DistFromOutside(const Double_t *point, const Double_t *dir, Int_t iact, Double_t step, Double_t *safe) const
493 {
494  // compute safe radius
495  if (iact<3 && safe) {
496  *safe = Safety(point, kFALSE);
497  if (iact==0) return TGeoShape::Big();
498  if ((iact==1) && (*safe>step)) return TGeoShape::Big();
499  }
500  // Check if the bounding box is crossed within the requested distance
501  Double_t sdist = TGeoBBox::DistFromOutside(point,dir, fDX, fDY, fDZ, fOrigin, step);
502  if (sdist>=step) return TGeoShape::Big();
503  // compute distance to Z planes
504  return TGeoCone::DistFromOutsideS(point, dir, fDz, fRmin1, fRmax1, fRmin2, fRmax2);
505 }
506 
507 ////////////////////////////////////////////////////////////////////////////////
508 /// Static method to compute distance to a conical surface with :
509 /// - r1, z1 - radius and Z position of lower base
510 /// - r2, z2 - radius and Z position of upper base
511 
512 void TGeoCone::DistToCone(const Double_t *point, const Double_t *dir, Double_t dz, Double_t r1, Double_t r2,
513  Double_t &b, Double_t &delta)
514 {
515  delta = -1.;
516  if (dz<0) return;
517  Double_t ro0 = 0.5*(r1+r2);
518  Double_t tz = 0.5*(r2-r1)/dz;
519  Double_t rsq = point[0]*point[0] + point[1]*point[1];
520  Double_t rc = ro0 + point[2]*tz;
521 
522  Double_t a = dir[0]*dir[0] + dir[1]*dir[1] - tz*tz*dir[2]*dir[2];
523  b = point[0]*dir[0] + point[1]*dir[1] - tz*rc*dir[2];
524  Double_t c = rsq - rc*rc;
525 
526  if (TMath::Abs(a)<TGeoShape::Tolerance()) {
527  if (TMath::Abs(b)<TGeoShape::Tolerance()) return;
528  b = 0.5*c/b;
529  delta = 0.;
530  return;
531  }
532  a = 1./a;
533  b *= a;
534  c *= a;
535  delta = b*b - c;
536  if (delta>0) {
537  delta = TMath::Sqrt(delta);
538  } else {
539  delta = -1.;
540  }
541 }
542 
543 ////////////////////////////////////////////////////////////////////////////////
544 /// compute closest distance from point px,py to each corner
545 
546 Int_t TGeoCone::DistancetoPrimitive(Int_t px, Int_t py)
547 {
548  Int_t n = gGeoManager->GetNsegments();
549  const Int_t numPoints = 4*n;
550  return ShapeDistancetoPrimitive(numPoints, px, py);
551 }
552 
553 ////////////////////////////////////////////////////////////////////////////////
554 /// Divide this cone shape belonging to volume "voldiv" into ndiv volumes
555 /// called divname, from start position with the given step. Returns pointer
556 /// to created division cell volume in case of Z divisions. For Z division
557 /// creates all volumes with different shapes and returns pointer to volume that
558 /// was divided. In case a wrong division axis is supplied, returns pointer to
559 /// volume that was divided.
560 
561 TGeoVolume *TGeoCone::Divide(TGeoVolume *voldiv, const char *divname, Int_t iaxis, Int_t ndiv,
562  Double_t start, Double_t step)
563 {
564  TGeoShape *shape; //--- shape to be created
565  TGeoVolume *vol; //--- division volume to be created
566  TGeoVolumeMulti *vmulti; //--- generic divided volume
567  TGeoPatternFinder *finder; //--- finder to be attached
568  TString opt = ""; //--- option to be attached
569  Int_t id;
570  Double_t end = start+ndiv*step;
571  switch (iaxis) {
572  case 1: //--- R division
573  Error("Divide","division of a cone on R not implemented");
574  return 0;
575  case 2: // --- Phi division
576  finder = new TGeoPatternCylPhi(voldiv, ndiv, start, end);
577  voldiv->SetFinder(finder);
578  finder->SetDivIndex(voldiv->GetNdaughters());
579  shape = new TGeoConeSeg(fDz, fRmin1, fRmax1, fRmin2, fRmax2, -step/2, step/2);
580  vol = new TGeoVolume(divname, shape, voldiv->GetMedium());
581  vmulti = gGeoManager->MakeVolumeMulti(divname, voldiv->GetMedium());
582  vmulti->AddVolume(vol);
583  opt = "Phi";
584  for (id=0; id<ndiv; id++) {
585  voldiv->AddNodeOffset(vol, id, start+id*step+step/2, opt.Data());
586  ((TGeoNodeOffset*)voldiv->GetNodes()->At(voldiv->GetNdaughters()-1))->SetFinder(finder);
587  }
588  return vmulti;
589  case 3: //--- Z division
590  vmulti = gGeoManager->MakeVolumeMulti(divname, voldiv->GetMedium());
591  finder = new TGeoPatternZ(voldiv, ndiv, start, end);
592  voldiv->SetFinder(finder);
593  finder->SetDivIndex(voldiv->GetNdaughters());
594  for (id=0; id<ndiv; id++) {
595  Double_t z1 = start+id*step;
596  Double_t z2 = start+(id+1)*step;
597  Double_t rmin1n = 0.5*(fRmin1*(fDz-z1)+fRmin2*(fDz+z1))/fDz;
598  Double_t rmax1n = 0.5*(fRmax1*(fDz-z1)+fRmax2*(fDz+z1))/fDz;
599  Double_t rmin2n = 0.5*(fRmin1*(fDz-z2)+fRmin2*(fDz+z2))/fDz;
600  Double_t rmax2n = 0.5*(fRmax1*(fDz-z2)+fRmax2*(fDz+z2))/fDz;
601  shape = new TGeoCone(0.5*step,rmin1n, rmax1n, rmin2n, rmax2n);
602  vol = new TGeoVolume(divname, shape, voldiv->GetMedium());
603  vmulti->AddVolume(vol);
604  opt = "Z";
605  voldiv->AddNodeOffset(vol, id, start+id*step+step/2, opt.Data());
606  ((TGeoNodeOffset*)voldiv->GetNodes()->At(voldiv->GetNdaughters()-1))->SetFinder(finder);
607  }
608  return vmulti;
609  default:
610  Error("Divide", "Wrong axis type for division");
611  return 0;
612  }
613 }
614 
615 ////////////////////////////////////////////////////////////////////////////////
616 /// Returns name of axis IAXIS.
617 
618 const char *TGeoCone::GetAxisName(Int_t iaxis) const
619 {
620  switch (iaxis) {
621  case 1:
622  return "R";
623  case 2:
624  return "PHI";
625  case 3:
626  return "Z";
627  default:
628  return "undefined";
629  }
630 }
631 
632 ////////////////////////////////////////////////////////////////////////////////
633 /// Get range of shape for a given axis.
634 
635 Double_t TGeoCone::GetAxisRange(Int_t iaxis, Double_t &xlo, Double_t &xhi) const
636 {
637  xlo = 0;
638  xhi = 0;
639  Double_t dx = 0;
640  switch (iaxis) {
641  case 2:
642  xlo = 0.;
643  xhi = 360.;
644  return 360.;
645  case 3:
646  xlo = -fDz;
647  xhi = fDz;
648  dx = xhi-xlo;
649  return dx;
650  }
651  return dx;
652 }
653 
654 ////////////////////////////////////////////////////////////////////////////////
655 /// Fill vector param[4] with the bounding cylinder parameters. The order
656 /// is the following : Rmin, Rmax, Phi1, Phi2, dZ
657 
658 void TGeoCone::GetBoundingCylinder(Double_t *param) const
659 {
660  param[0] = TMath::Min(fRmin1, fRmin2); // Rmin
661  param[0] *= param[0];
662  param[1] = TMath::Max(fRmax1, fRmax2); // Rmax
663  param[1] *= param[1];
664  param[2] = 0.; // Phi1
665  param[3] = 360.; // Phi2
666 }
667 
668 ////////////////////////////////////////////////////////////////////////////////
669 /// in case shape has some negative parameters, these has to be computed
670 /// in order to fit the mother
671 
672 TGeoShape *TGeoCone::GetMakeRuntimeShape(TGeoShape *mother, TGeoMatrix * /*mat*/) const
673 {
674  if (!TestShapeBit(kGeoRunTimeShape)) return 0;
675  if (!mother->TestShapeBit(kGeoCone)) {
676  Error("GetMakeRuntimeShape", "invalid mother");
677  return 0;
678  }
679  Double_t rmin1, rmax1, rmin2, rmax2, dz;
680  rmin1 = fRmin1;
681  rmax1 = fRmax1;
682  rmin2 = fRmin2;
683  rmax2 = fRmax2;
684  dz = fDz;
685  if (fDz<0) dz=((TGeoCone*)mother)->GetDz();
686  if (fRmin1<0)
687  rmin1 = ((TGeoCone*)mother)->GetRmin1();
688  if (fRmax1<0)
689  rmax1 = ((TGeoCone*)mother)->GetRmax1();
690  if (fRmin2<0)
691  rmin2 = ((TGeoCone*)mother)->GetRmin2();
692  if (fRmax2<0)
693  rmax2 = ((TGeoCone*)mother)->GetRmax2();
694 
695  return (new TGeoCone(GetName(), dz, rmin1, rmax1, rmin2, rmax2));
696 }
697 
698 ////////////////////////////////////////////////////////////////////////////////
699 /// Fills array with n random points located on the line segments of the shape mesh.
700 /// The output array must be provided with a length of minimum 3*npoints. Returns
701 /// true if operation is implemented.
702 
703 Bool_t TGeoCone::GetPointsOnSegments(Int_t npoints, Double_t *array) const
704 {
705  if (npoints > (npoints/2)*2) {
706  Error("GetPointsOnSegments","Npoints must be even number");
707  return kFALSE;
708  }
709  Bool_t hasrmin = (fRmin1>0 || fRmin2>0)?kTRUE:kFALSE;
710  Int_t nc = 0;
711  if (hasrmin) nc = (Int_t)TMath::Sqrt(0.5*npoints);
712  else nc = (Int_t)TMath::Sqrt(1.*npoints);
713  Double_t dphi = TMath::TwoPi()/nc;
714  Double_t phi = 0;
715  Int_t ntop = 0;
716  if (hasrmin) ntop = npoints/2 - nc*(nc-1);
717  else ntop = npoints - nc*(nc-1);
718  Double_t dz = 2*fDz/(nc-1);
719  Double_t z = 0;
720  Int_t icrt = 0;
721  Int_t nphi = nc;
722  Double_t rmin = 0.;
723  Double_t rmax = 0.;
724  // loop z sections
725  for (Int_t i=0; i<nc; i++) {
726  if (i == (nc-1)) nphi = ntop;
727  z = -fDz + i*dz;
728  if (hasrmin) rmin = 0.5*(fRmin1+fRmin2) + 0.5*(fRmin2-fRmin1)*z/fDz;
729  rmax = 0.5*(fRmax1+fRmax2) + 0.5*(fRmax2-fRmax1)*z/fDz;
730  // loop points on circle sections
731  for (Int_t j=0; j<nphi; j++) {
732  phi = j*dphi;
733  if (hasrmin) {
734  array[icrt++] = rmin * TMath::Cos(phi);
735  array[icrt++] = rmin * TMath::Sin(phi);
736  array[icrt++] = z;
737  }
738  array[icrt++] = rmax * TMath::Cos(phi);
739  array[icrt++] = rmax * TMath::Sin(phi);
740  array[icrt++] = z;
741  }
742  }
743  return kTRUE;
744 }
745 
746 
747 ////////////////////////////////////////////////////////////////////////////////
748 /// print shape parameters
749 
750 void TGeoCone::InspectShape() const
751 {
752  printf("*** Shape %s TGeoCone ***\n", GetName());
753  printf(" dz =: %11.5f\n", fDz);
754  printf(" Rmin1 = %11.5f\n", fRmin1);
755  printf(" Rmax1 = %11.5f\n", fRmax1);
756  printf(" Rmin2 = %11.5f\n", fRmin2);
757  printf(" Rmax2 = %11.5f\n", fRmax2);
758  printf(" Bounding box:\n");
759  TGeoBBox::InspectShape();
760 }
761 
762 ////////////////////////////////////////////////////////////////////////////////
763 /// Creates a TBuffer3D describing *this* shape.
764 /// Coordinates are in local reference frame.
765 
766 TBuffer3D *TGeoCone::MakeBuffer3D() const
767 {
768  Int_t n = gGeoManager->GetNsegments();
769  Int_t nbPnts = 4*n;
770  Int_t nbSegs = 8*n;
771  Int_t nbPols = 4*n;
772  TBuffer3D* buff = new TBuffer3D(TBuffer3DTypes::kGeneric,
773  nbPnts, 3*nbPnts,
774  nbSegs, 3*nbSegs,
775  nbPols, 6*nbPols);
776 
777  if (buff)
778  {
779  SetPoints(buff->fPnts);
780  SetSegsAndPols(*buff);
781  }
782  return buff;
783 }
784 
785 ////////////////////////////////////////////////////////////////////////////////
786 /// Fill TBuffer3D structure for segments and polygons.
787 
788 void TGeoCone::SetSegsAndPols(TBuffer3D &buffer) const
789 {
790  Int_t i,j;
791  Int_t n = gGeoManager->GetNsegments();
792  Int_t c = GetBasicColor();
793 
794  for (i = 0; i < 4; i++) {
795  for (j = 0; j < n; j++) {
796  buffer.fSegs[(i*n+j)*3 ] = c;
797  buffer.fSegs[(i*n+j)*3+1] = i*n+j;
798  buffer.fSegs[(i*n+j)*3+2] = i*n+j+1;
799  }
800  buffer.fSegs[(i*n+j-1)*3+2] = i*n;
801  }
802  for (i = 4; i < 6; i++) {
803  for (j = 0; j < n; j++) {
804  buffer.fSegs[(i*n+j)*3 ] = c+1;
805  buffer.fSegs[(i*n+j)*3+1] = (i-4)*n+j;
806  buffer.fSegs[(i*n+j)*3+2] = (i-2)*n+j;
807  }
808  }
809  for (i = 6; i < 8; i++) {
810  for (j = 0; j < n; j++) {
811  buffer.fSegs[(i*n+j)*3 ] = c;
812  buffer.fSegs[(i*n+j)*3+1] = 2*(i-6)*n+j;
813  buffer.fSegs[(i*n+j)*3+2] = (2*(i-6)+1)*n+j;
814  }
815  }
816 
817  Int_t indx = 0;
818  i=0;
819  for (j = 0; j < n; j++) {
820  indx = 6*(i*n+j);
821  buffer.fPols[indx ] = c;
822  buffer.fPols[indx+1] = 4;
823  buffer.fPols[indx+5] = i*n+j;
824  buffer.fPols[indx+4] = (4+i)*n+j;
825  buffer.fPols[indx+3] = (2+i)*n+j;
826  buffer.fPols[indx+2] = (4+i)*n+j+1;
827  }
828  buffer.fPols[indx+2] = (4+i)*n;
829  i=1;
830  for (j = 0; j < n; j++) {
831  indx = 6*(i*n+j);
832  buffer.fPols[indx ] = c;
833  buffer.fPols[indx+1] = 4;
834  buffer.fPols[indx+2] = i*n+j;
835  buffer.fPols[indx+3] = (4+i)*n+j;
836  buffer.fPols[indx+4] = (2+i)*n+j;
837  buffer.fPols[indx+5] = (4+i)*n+j+1;
838  }
839  buffer.fPols[indx+5] = (4+i)*n;
840  i=2;
841  for (j = 0; j < n; j++) {
842  indx = 6*(i*n+j);
843  buffer.fPols[indx ] = c+i;
844  buffer.fPols[indx+1] = 4;
845  buffer.fPols[indx+2] = (i-2)*2*n+j;
846  buffer.fPols[indx+3] = (4+i)*n+j;
847  buffer.fPols[indx+4] = ((i-2)*2+1)*n+j;
848  buffer.fPols[indx+5] = (4+i)*n+j+1;
849  }
850  buffer.fPols[indx+5] = (4+i)*n;
851  i=3;
852  for (j = 0; j < n; j++) {
853  indx = 6*(i*n+j);
854  buffer.fPols[indx ] = c+i;
855  buffer.fPols[indx+1] = 4;
856  buffer.fPols[indx+5] = (i-2)*2*n+j;
857  buffer.fPols[indx+4] = (4+i)*n+j;
858  buffer.fPols[indx+3] = ((i-2)*2+1)*n+j;
859  buffer.fPols[indx+2] = (4+i)*n+j+1;
860  }
861  buffer.fPols[indx+2] = (4+i)*n;
862 }
863 
864 ////////////////////////////////////////////////////////////////////////////////
865 /// computes the closest distance from given point to this shape, according
866 /// to option. The matching point on the shape is stored in spoint.
867 
868 Double_t TGeoCone::Safety(const Double_t *point, Bool_t in) const
869 {
870  Double_t saf[4];
871  Double_t r=TMath::Sqrt(point[0]*point[0]+point[1]*point[1]);
872  saf[0] = TGeoShape::SafetySeg(r,point[2], fRmin1, -fDz, fRmax1, -fDz, !in);
873  saf[1] = TGeoShape::SafetySeg(r,point[2], fRmax2, fDz, fRmin2, fDz, !in);
874  saf[2] = TGeoShape::SafetySeg(r,point[2], fRmin2, fDz, fRmin1, -fDz, !in);
875  saf[3] = TGeoShape::SafetySeg(r,point[2], fRmax1, -fDz, fRmax2, fDz, !in);
876  Double_t safety = saf[TMath::LocMin(4,saf)];
877  if (safety>1.E20) safety = 0.;
878  return safety;
879 }
880 
881 ////////////////////////////////////////////////////////////////////////////////
882 /// computes the closest distance from given point to this shape, according
883 /// to option. The matching point on the shape is stored in spoint.
884 
885 Double_t TGeoCone::SafetyS(const Double_t *point, Bool_t in, Double_t dz, Double_t rmin1, Double_t rmax1,
886  Double_t rmin2, Double_t rmax2, Int_t skipz)
887 {
888  Double_t saf[4];
889  Double_t r=TMath::Sqrt(point[0]*point[0]+point[1]*point[1]);
890 // Double_t rin = tg1*point[2]+ro1;
891 // Double_t rout = tg2*point[2]+ro2;
892  switch (skipz) {
893  case 1: // skip lower Z plane
894  saf[0] = TGeoShape::Big();
895  saf[1] = TGeoShape::SafetySeg(r,point[2], rmax2, dz, rmin2, dz, !in);
896  break;
897  case 2: // skip upper Z plane
898  saf[0] = TGeoShape::SafetySeg(r,point[2], rmin1, -dz, rmax1, -dz, !in);
899  saf[1] = TGeoShape::Big();
900  break;
901  case 3: // skip both
902  saf[0] = saf[1] = TGeoShape::Big();
903  break;
904  default:
905  saf[0] = TGeoShape::SafetySeg(r,point[2], rmin1, -dz, rmax1, -dz, !in);
906  saf[1] = TGeoShape::SafetySeg(r,point[2], rmax2, dz, rmin2, dz, !in);
907  }
908  // Safety to inner part
909  if (rmin1>0 || rmin2>0)
910  saf[2] = TGeoShape::SafetySeg(r,point[2], rmin2, dz, rmin1, -dz, !in);
911  else
912  saf[2] = TGeoShape::Big();
913  saf[3] = TGeoShape::SafetySeg(r,point[2], rmax1, -dz, rmax2, dz, !in);
914  return saf[TMath::LocMin(4,saf)];
915 }
916 
917 ////////////////////////////////////////////////////////////////////////////////
918 /// Save a primitive as a C++ statement(s) on output stream "out".
919 
920 void TGeoCone::SavePrimitive(std::ostream &out, Option_t * /*option*/ /*= ""*/)
921 {
922  if (TObject::TestBit(kGeoSavePrimitive)) return;
923  out << " // Shape: " << GetName() << " type: " << ClassName() << std::endl;
924  out << " dz = " << fDz << ";" << std::endl;
925  out << " rmin1 = " << fRmin1 << ";" << std::endl;
926  out << " rmax1 = " << fRmax1 << ";" << std::endl;
927  out << " rmin2 = " << fRmin2 << ";" << std::endl;
928  out << " rmax2 = " << fRmax2 << ";" << std::endl;
929  out << " TGeoShape *" << GetPointerName() << " = new TGeoCone(\"" << GetName() << "\", dz,rmin1,rmax1,rmin2,rmax2);" << std::endl;
930  TObject::SetBit(TGeoShape::kGeoSavePrimitive);
931 }
932 
933 ////////////////////////////////////////////////////////////////////////////////
934 /// Set cone dimensions.
935 
936 void TGeoCone::SetConeDimensions(Double_t dz, Double_t rmin1, Double_t rmax1,
937  Double_t rmin2, Double_t rmax2)
938 {
939  if (rmin1>=0) {
940  if (rmax1>0) {
941  if (rmin1<=rmax1) {
942  // normal rmin/rmax
943  fRmin1 = rmin1;
944  fRmax1 = rmax1;
945  } else {
946  fRmin1 = rmax1;
947  fRmax1 = rmin1;
948  Warning("SetConeDimensions", "rmin1>rmax1 Switch rmin1<->rmax1");
949  SetShapeBit(TGeoShape::kGeoBad);
950  }
951  } else {
952  // run-time
953  fRmin1 = rmin1;
954  fRmax1 = rmax1;
955  }
956  } else {
957  // run-time
958  fRmin1 = rmin1;
959  fRmax1 = rmax1;
960  }
961  if (rmin2>=0) {
962  if (rmax2>0) {
963  if (rmin2<=rmax2) {
964  // normal rmin/rmax
965  fRmin2 = rmin2;
966  fRmax2 = rmax2;
967  } else {
968  fRmin2 = rmax2;
969  fRmax2 = rmin2;
970  Warning("SetConeDimensions", "rmin2>rmax2 Switch rmin2<->rmax2");
971  SetShapeBit(TGeoShape::kGeoBad);
972  }
973  } else {
974  // run-time
975  fRmin2 = rmin2;
976  fRmax2 = rmax2;
977  }
978  } else {
979  // run-time
980  fRmin2 = rmin2;
981  fRmax2 = rmax2;
982  }
983 
984  fDz = dz;
985 }
986 
987 ////////////////////////////////////////////////////////////////////////////////
988 /// Set cone dimensions from an array.
989 
990 void TGeoCone::SetDimensions(Double_t *param)
991 {
992  Double_t dz = param[0];
993  Double_t rmin1 = param[1];
994  Double_t rmax1 = param[2];
995  Double_t rmin2 = param[3];
996  Double_t rmax2 = param[4];
997  SetConeDimensions(dz, rmin1, rmax1, rmin2, rmax2);
998 }
999 
1000 ////////////////////////////////////////////////////////////////////////////////
1001 /// Create cone mesh points.
1002 
1003 void TGeoCone::SetPoints(Double_t *points) const
1004 {
1005  Double_t dz, phi, dphi;
1006  Int_t j, n;
1007 
1008  n = gGeoManager->GetNsegments();
1009  dphi = 360./n;
1010  dz = fDz;
1011  Int_t indx = 0;
1012 
1013  if (points) {
1014  for (j = 0; j < n; j++) {
1015  phi = j*dphi*TMath::DegToRad();
1016  points[indx++] = fRmin1 * TMath::Cos(phi);
1017  points[indx++] = fRmin1 * TMath::Sin(phi);
1018  points[indx++] = -dz;
1019  }
1020 
1021  for (j = 0; j < n; j++) {
1022  phi = j*dphi*TMath::DegToRad();
1023  points[indx++] = fRmax1 * TMath::Cos(phi);
1024  points[indx++] = fRmax1 * TMath::Sin(phi);
1025  points[indx++] = -dz;
1026  }
1027 
1028  for (j = 0; j < n; j++) {
1029  phi = j*dphi*TMath::DegToRad();
1030  points[indx++] = fRmin2 * TMath::Cos(phi);
1031  points[indx++] = fRmin2 * TMath::Sin(phi);
1032  points[indx++] = dz;
1033  }
1034 
1035  for (j = 0; j < n; j++) {
1036  phi = j*dphi*TMath::DegToRad();
1037  points[indx++] = fRmax2 * TMath::Cos(phi);
1038  points[indx++] = fRmax2 * TMath::Sin(phi);
1039  points[indx++] = dz;
1040  }
1041  }
1042 }
1043 
1044 ////////////////////////////////////////////////////////////////////////////////
1045 /// Create cone mesh points.
1046 
1047 void TGeoCone::SetPoints(Float_t *points) const
1048 {
1049  Double_t dz, phi, dphi;
1050  Int_t j, n;
1051 
1052  n = gGeoManager->GetNsegments();
1053  dphi = 360./n;
1054  dz = fDz;
1055  Int_t indx = 0;
1056 
1057  if (points) {
1058  for (j = 0; j < n; j++) {
1059  phi = j*dphi*TMath::DegToRad();
1060  points[indx++] = fRmin1 * TMath::Cos(phi);
1061  points[indx++] = fRmin1 * TMath::Sin(phi);
1062  points[indx++] = -dz;
1063  }
1064 
1065  for (j = 0; j < n; j++) {
1066  phi = j*dphi*TMath::DegToRad();
1067  points[indx++] = fRmax1 * TMath::Cos(phi);
1068  points[indx++] = fRmax1 * TMath::Sin(phi);
1069  points[indx++] = -dz;
1070  }
1071 
1072  for (j = 0; j < n; j++) {
1073  phi = j*dphi*TMath::DegToRad();
1074  points[indx++] = fRmin2 * TMath::Cos(phi);
1075  points[indx++] = fRmin2 * TMath::Sin(phi);
1076  points[indx++] = dz;
1077  }
1078 
1079  for (j = 0; j < n; j++) {
1080  phi = j*dphi*TMath::DegToRad();
1081  points[indx++] = fRmax2 * TMath::Cos(phi);
1082  points[indx++] = fRmax2 * TMath::Sin(phi);
1083  points[indx++] = dz;
1084  }
1085  }
1086 }
1087 
1088 ////////////////////////////////////////////////////////////////////////////////
1089 /// Returns numbers of vertices, segments and polygons composing the shape mesh.
1090 
1091 void TGeoCone::GetMeshNumbers(Int_t &nvert, Int_t &nsegs, Int_t &npols) const
1092 {
1093  Int_t n = gGeoManager->GetNsegments();
1094  nvert = n*4;
1095  nsegs = n*8;
1096  npols = n*4;
1097 }
1098 
1099 ////////////////////////////////////////////////////////////////////////////////
1100 /// Return number of vertices of the mesh representation
1101 
1102 Int_t TGeoCone::GetNmeshVertices() const
1103 {
1104  Int_t n = gGeoManager->GetNsegments();
1105  Int_t numPoints = n*4;
1106  return numPoints;
1107 }
1108 
1109 ////////////////////////////////////////////////////////////////////////////////
1110 /// Fill size of this 3-D object
1111 
1112 void TGeoCone::Sizeof3D() const
1113 {
1114 }
1115 
1116 ////////////////////////////////////////////////////////////////////////////////
1117 /// Fills a static 3D buffer and returns a reference.
1118 
1119 const TBuffer3D & TGeoCone::GetBuffer3D(Int_t reqSections, Bool_t localFrame) const
1120 {
1121  static TBuffer3D buffer(TBuffer3DTypes::kGeneric);
1122 
1123  TGeoBBox::FillBuffer3D(buffer, reqSections, localFrame);
1124 
1125  if (reqSections & TBuffer3D::kRawSizes) {
1126  Int_t n = gGeoManager->GetNsegments();
1127  Int_t nbPnts = 4*n;
1128  Int_t nbSegs = 8*n;
1129  Int_t nbPols = 4*n;
1130  if (buffer.SetRawSizes(nbPnts, 3*nbPnts, nbSegs, 3*nbSegs, nbPols, 6*nbPols)) {
1131  buffer.SetSectionsValid(TBuffer3D::kRawSizes);
1132  }
1133  }
1134 
1135  // TODO: Can we push this as common down to TGeoShape?
1136  if ((reqSections & TBuffer3D::kRaw) && buffer.SectionsValid(TBuffer3D::kRawSizes)) {
1137  SetPoints(buffer.fPnts);
1138  if (!buffer.fLocalFrame) {
1139  TransformPoints(buffer.fPnts, buffer.NbPnts());
1140  }
1141 
1142  SetSegsAndPols(buffer);
1143  buffer.SetSectionsValid(TBuffer3D::kRaw);
1144  }
1145 
1146  return buffer;
1147 }
1148 
1149 ////////////////////////////////////////////////////////////////////////////////
1150 /// Check the inside status for each of the points in the array.
1151 /// Input: Array of point coordinates + vector size
1152 /// Output: Array of Booleans for the inside of each point
1153 
1154 void TGeoCone::Contains_v(const Double_t *points, Bool_t *inside, Int_t vecsize) const
1155 {
1156  for (Int_t i=0; i<vecsize; i++) inside[i] = Contains(&points[3*i]);
1157 }
1158 
1159 ////////////////////////////////////////////////////////////////////////////////
1160 /// Compute the normal for an array o points so that norm.dot.dir is positive
1161 /// Input: Arrays of point coordinates and directions + vector size
1162 /// Output: Array of normal directions
1163 
1164 void TGeoCone::ComputeNormal_v(const Double_t *points, const Double_t *dirs, Double_t *norms, Int_t vecsize)
1165 {
1166  for (Int_t i=0; i<vecsize; i++) ComputeNormal(&points[3*i], &dirs[3*i], &norms[3*i]);
1167 }
1168 
1169 ////////////////////////////////////////////////////////////////////////////////
1170 /// Compute distance from array of input points having directions specified by dirs. Store output in dists
1171 
1172 void TGeoCone::DistFromInside_v(const Double_t *points, const Double_t *dirs, Double_t *dists, Int_t vecsize, Double_t* step) const
1173 {
1174  for (Int_t i=0; i<vecsize; i++) dists[i] = DistFromInside(&points[3*i], &dirs[3*i], 3, step[i]);
1175 }
1176 
1177 ////////////////////////////////////////////////////////////////////////////////
1178 /// Compute distance from array of input points having directions specified by dirs. Store output in dists
1179 
1180 void TGeoCone::DistFromOutside_v(const Double_t *points, const Double_t *dirs, Double_t *dists, Int_t vecsize, Double_t* step) const
1181 {
1182  for (Int_t i=0; i<vecsize; i++) dists[i] = DistFromOutside(&points[3*i], &dirs[3*i], 3, step[i]);
1183 }
1184 
1185 ////////////////////////////////////////////////////////////////////////////////
1186 /// Compute safe distance from each of the points in the input array.
1187 /// Input: Array of point coordinates, array of statuses for these points, size of the arrays
1188 /// Output: Safety values
1189 
1190 void TGeoCone::Safety_v(const Double_t *points, const Bool_t *inside, Double_t *safe, Int_t vecsize) const
1191 {
1192  for (Int_t i=0; i<vecsize; i++) safe[i] = Safety(&points[3*i], inside[i]);
1193 }
1194 
1195 ClassImp(TGeoConeSeg);
1196 
1197 ////////////////////////////////////////////////////////////////////////////////
1198 /// Default constructor
1199 
1200 TGeoConeSeg::TGeoConeSeg()
1201  :TGeoCone(),
1202  fPhi1(0.), fPhi2(0.), fS1(0.), fC1(0.), fS2(0.), fC2(0.), fSm(0.), fCm(0.), fCdfi(0.)
1203 {
1204  SetShapeBit(TGeoShape::kGeoConeSeg);
1205  fPhi1 = fPhi2 = 0.0;
1206 }
1207 
1208 ////////////////////////////////////////////////////////////////////////////////
1209 /// Default constructor specifying minimum and maximum radius
1210 
1211 TGeoConeSeg::TGeoConeSeg(Double_t dz, Double_t rmin1, Double_t rmax1,
1212  Double_t rmin2, Double_t rmax2, Double_t phi1, Double_t phi2)
1213  :TGeoCone(dz, rmin1, rmax1, rmin2, rmax2),
1214  fPhi1(0.), fPhi2(0.), fS1(0.), fC1(0.), fS2(0.), fC2(0.), fSm(0.), fCm(0.), fCdfi(0.)
1215 
1216 {
1217  SetShapeBit(TGeoShape::kGeoConeSeg);
1218  SetConsDimensions(dz, rmin1, rmax1, rmin2, rmax2, phi1, phi2);
1219  ComputeBBox();
1220 }
1221 
1222 ////////////////////////////////////////////////////////////////////////////////
1223 /// Default constructor specifying minimum and maximum radius
1224 
1225 TGeoConeSeg::TGeoConeSeg(const char *name, Double_t dz, Double_t rmin1, Double_t rmax1,
1226  Double_t rmin2, Double_t rmax2, Double_t phi1, Double_t phi2)
1227  :TGeoCone(name, dz, rmin1, rmax1, rmin2, rmax2),
1228  fPhi1(0.), fPhi2(0.), fS1(0.), fC1(0.), fS2(0.), fC2(0.), fSm(0.), fCm(0.), fCdfi(0.)
1229 {
1230  SetShapeBit(TGeoShape::kGeoConeSeg);
1231  SetConsDimensions(dz, rmin1, rmax1, rmin2, rmax2, phi1, phi2);
1232  ComputeBBox();
1233 }
1234 
1235 ////////////////////////////////////////////////////////////////////////////////
1236 /// Default constructor specifying minimum and maximum radius
1237 /// - param[0] = dz
1238 /// - param[1] = Rmin1
1239 /// - param[2] = Rmax1
1240 /// - param[3] = Rmin2
1241 /// - param[4] = Rmax2
1242 /// - param[5] = phi1
1243 /// - param[6] = phi2
1244 
1245 TGeoConeSeg::TGeoConeSeg(Double_t *param)
1246  :TGeoCone(0,0,0,0,0),
1247  fPhi1(0.), fPhi2(0.), fS1(0.), fC1(0.), fS2(0.), fC2(0.), fSm(0.), fCm(0.), fCdfi(0.)
1248 {
1249  SetShapeBit(TGeoShape::kGeoConeSeg);
1250  SetDimensions(param);
1251  ComputeBBox();
1252 }
1253 
1254 ////////////////////////////////////////////////////////////////////////////////
1255 /// destructor
1256 
1257 TGeoConeSeg::~TGeoConeSeg()
1258 {
1259 }
1260 
1261 ////////////////////////////////////////////////////////////////////////////////
1262 /// Function called after streaming an object of this class.
1263 
1264 void TGeoConeSeg::AfterStreamer()
1265 {
1266  InitTrigonometry();
1267 }
1268 
1269 ////////////////////////////////////////////////////////////////////////////////
1270 /// Init frequently used trigonometric values
1271 
1272 void TGeoConeSeg::InitTrigonometry()
1273 {
1274  Double_t phi1 = fPhi1*TMath::DegToRad();
1275  Double_t phi2 = fPhi2*TMath::DegToRad();
1276  fC1 = TMath::Cos(phi1);
1277  fS1 = TMath::Sin(phi1);
1278  fC2 = TMath::Cos(phi2);
1279  fS2 = TMath::Sin(phi2);
1280  Double_t fio = 0.5*(phi1+phi2);
1281  fCm = TMath::Cos(fio);
1282  fSm = TMath::Sin(fio);
1283  Double_t dfi = 0.5*(phi2-phi1);
1284  fCdfi = TMath::Cos(dfi);
1285 }
1286 
1287 ////////////////////////////////////////////////////////////////////////////////
1288 /// Computes capacity of the shape in [length^3]
1289 
1290 Double_t TGeoConeSeg::Capacity() const
1291 {
1292  return TGeoConeSeg::Capacity(fDz, fRmin1, fRmax1, fRmin2, fRmax2, fPhi1, fPhi2);
1293 }
1294 
1295 ////////////////////////////////////////////////////////////////////////////////
1296 /// Computes capacity of the shape in [length^3]
1297 
1298 Double_t TGeoConeSeg::Capacity(Double_t dz, Double_t rmin1, Double_t rmax1, Double_t rmin2, Double_t rmax2, Double_t phi1, Double_t phi2)
1299 {
1300  Double_t capacity = (TMath::Abs(phi2-phi1)*TMath::DegToRad()*dz/3.)*
1301  (rmax1*rmax1+rmax2*rmax2+rmax1*rmax2-
1302  rmin1*rmin1-rmin2*rmin2-rmin1*rmin2);
1303  return capacity;
1304 }
1305 
1306 ////////////////////////////////////////////////////////////////////////////////
1307 /// compute bounding box of the tube segment
1308 
1309 void TGeoConeSeg::ComputeBBox()
1310 {
1311  Double_t rmin, rmax;
1312  rmin = TMath::Min(fRmin1, fRmin2);
1313  rmax = TMath::Max(fRmax1, fRmax2);
1314 
1315  Double_t xc[4];
1316  Double_t yc[4];
1317  xc[0] = rmax*fC1;
1318  yc[0] = rmax*fS1;
1319  xc[1] = rmax*fC2;
1320  yc[1] = rmax*fS2;
1321  xc[2] = rmin*fC1;
1322  yc[2] = rmin*fS1;
1323  xc[3] = rmin*fC2;
1324  yc[3] = rmin*fS2;
1325 
1326  Double_t xmin = xc[TMath::LocMin(4, &xc[0])];
1327  Double_t xmax = xc[TMath::LocMax(4, &xc[0])];
1328  Double_t ymin = yc[TMath::LocMin(4, &yc[0])];
1329  Double_t ymax = yc[TMath::LocMax(4, &yc[0])];
1330 
1331  Double_t dp = fPhi2-fPhi1;
1332  Double_t ddp = -fPhi1;
1333  if (ddp<0) ddp+= 360;
1334  if (ddp<=dp) xmax = rmax;
1335  ddp = 90-fPhi1;
1336  if (ddp<0) ddp+= 360;
1337  if (ddp<=dp) ymax = rmax;
1338  ddp = 180-fPhi1;
1339  if (ddp<0) ddp+= 360;
1340  if (ddp<=dp) xmin = -rmax;
1341  ddp = 270-fPhi1;
1342  if (ddp<0) ddp+= 360;
1343  if (ddp<=dp) ymin = -rmax;
1344  fOrigin[0] = (xmax+xmin)/2;
1345  fOrigin[1] = (ymax+ymin)/2;
1346  fOrigin[2] = 0;
1347  fDX = (xmax-xmin)/2;
1348  fDY = (ymax-ymin)/2;
1349  fDZ = fDz;
1350 }
1351 
1352 ////////////////////////////////////////////////////////////////////////////////
1353 /// Compute normal to closest surface from POINT.
1354 
1355 void TGeoConeSeg::ComputeNormal(const Double_t *point, const Double_t *dir, Double_t *norm)
1356 {
1357  Double_t saf[3];
1358  Double_t ro1 = 0.5*(fRmin1+fRmin2);
1359  Double_t tg1 = 0.5*(fRmin2-fRmin1)/fDz;
1360  Double_t cr1 = 1./TMath::Sqrt(1.+tg1*tg1);
1361  Double_t ro2 = 0.5*(fRmax1+fRmax2);
1362  Double_t tg2 = 0.5*(fRmax2-fRmax1)/fDz;
1363  Double_t cr2 = 1./TMath::Sqrt(1.+tg2*tg2);
1364 
1365  Double_t r=TMath::Sqrt(point[0]*point[0]+point[1]*point[1]);
1366  Double_t rin = tg1*point[2]+ro1;
1367  Double_t rout = tg2*point[2]+ro2;
1368  saf[0] = TMath::Abs(fDz-TMath::Abs(point[2]));
1369  saf[1] = (ro1>0)?(TMath::Abs((r-rin)*cr1)):TGeoShape::Big();
1370  saf[2] = TMath::Abs((rout-r)*cr2);
1371  Int_t i = TMath::LocMin(3,saf);
1372  if (((fPhi2-fPhi1)<360.) && TGeoShape::IsCloseToPhi(saf[i], point,fC1,fS1,fC2,fS2)) {
1373  TGeoShape::NormalPhi(point,dir,norm,fC1,fS1,fC2,fS2);
1374  return;
1375  }
1376  if (i==0) {
1377  norm[0] = norm[1] = 0.;
1378  norm[2] = TMath::Sign(1.,dir[2]);
1379  return;
1380  }
1381 
1382  Double_t phi = TMath::ATan2(point[1],point[0]);
1383  Double_t cphi = TMath::Cos(phi);
1384  Double_t sphi = TMath::Sin(phi);
1385 
1386  if (i==1) {
1387  norm[0] = cr1*cphi;
1388  norm[1] = cr1*sphi;
1389  norm[2] = -tg1*cr1;
1390  } else {
1391  norm[0] = cr2*cphi;
1392  norm[1] = cr2*sphi;
1393  norm[2] = -tg2*cr2;
1394  }
1395 
1396  if (norm[0]*dir[0]+norm[1]*dir[1]+norm[2]*dir[2]<0) {
1397  norm[0] = -norm[0];
1398  norm[1] = -norm[1];
1399  norm[2] = -norm[2];
1400  }
1401 }
1402 
1403 ////////////////////////////////////////////////////////////////////////////////
1404 /// Compute normal to closest surface from POINT.
1405 
1406 void TGeoConeSeg::ComputeNormalS(const Double_t *point, const Double_t *dir, Double_t *norm,
1407  Double_t dz, Double_t rmin1, Double_t rmax1, Double_t rmin2, Double_t rmax2,
1408  Double_t c1, Double_t s1, Double_t c2, Double_t s2)
1409 {
1410  Double_t saf[2];
1411  Double_t ro1 = 0.5*(rmin1+rmin2);
1412  Double_t tg1 = 0.5*(rmin2-rmin1)/dz;
1413  Double_t cr1 = 1./TMath::Sqrt(1.+tg1*tg1);
1414  Double_t ro2 = 0.5*(rmax1+rmax2);
1415  Double_t tg2 = 0.5*(rmax2-rmax1)/dz;
1416  Double_t cr2 = 1./TMath::Sqrt(1.+tg2*tg2);
1417 
1418  Double_t r=TMath::Sqrt(point[0]*point[0]+point[1]*point[1]);
1419  Double_t rin = tg1*point[2]+ro1;
1420  Double_t rout = tg2*point[2]+ro2;
1421  saf[0] = (ro1>0)?(TMath::Abs((r-rin)*cr1)):TGeoShape::Big();
1422  saf[1] = TMath::Abs((rout-r)*cr2);
1423  Int_t i = TMath::LocMin(2,saf);
1424  if (TGeoShape::IsCloseToPhi(saf[i], point,c1,s1,c2,s2)) {
1425  TGeoShape::NormalPhi(point,dir,norm,c1,s1,c2,s2);
1426  return;
1427  }
1428 
1429  Double_t phi = TMath::ATan2(point[1],point[0]);
1430  Double_t cphi = TMath::Cos(phi);
1431  Double_t sphi = TMath::Sin(phi);
1432 
1433  if (i==0) {
1434  norm[0] = cr1*cphi;
1435  norm[1] = cr1*sphi;
1436  norm[2] = -tg1*cr1;
1437  } else {
1438  norm[0] = cr2*cphi;
1439  norm[1] = cr2*sphi;
1440  norm[2] = -tg2*cr2;
1441  }
1442 
1443  if (norm[0]*dir[0]+norm[1]*dir[1]+norm[2]*dir[2]<0) {
1444  norm[0] = -norm[0];
1445  norm[1] = -norm[1];
1446  norm[2] = -norm[2];
1447  }
1448 }
1449 
1450 ////////////////////////////////////////////////////////////////////////////////
1451 /// test if point is inside this sphere
1452 
1453 Bool_t TGeoConeSeg::Contains(const Double_t *point) const
1454 {
1455  if (!TGeoCone::Contains(point)) return kFALSE;
1456  Double_t dphi = fPhi2 - fPhi1;
1457  if (dphi >= 360.) return kTRUE;
1458  Double_t phi = TMath::ATan2(point[1], point[0]) * TMath::RadToDeg();
1459  if (phi < 0 ) phi+=360.;
1460  Double_t ddp = phi-fPhi1;
1461  if (ddp < 0) ddp+=360.;
1462 // if (ddp > 360) ddp-=360;
1463  if (ddp > dphi) return kFALSE;
1464  return kTRUE;
1465 }
1466 
1467 ////////////////////////////////////////////////////////////////////////////////
1468 /// Static method to compute distance to a conical surface with :
1469 /// - r1, z1 - radius and Z position of lower base
1470 /// - r2, z2 - radius and Z position of upper base
1471 /// - phi1, phi2 - phi limits
1472 
1473 Double_t TGeoConeSeg::DistToCons(const Double_t *point, const Double_t *dir, Double_t r1, Double_t z1, Double_t r2, Double_t z2, Double_t phi1, Double_t phi2)
1474 {
1475  Double_t dz = z2-z1;
1476  if (dz<=0) {
1477  return TGeoShape::Big();
1478  }
1479 
1480  Double_t dphi = phi2 - phi1;
1481  Bool_t hasphi = kTRUE;
1482  if (dphi >= 360.) hasphi=kFALSE;
1483  if (dphi < 0) dphi+=360.;
1484 // printf("phi1=%f phi2=%f dphi=%f\n", phi1, phi2, dphi);
1485 
1486  Double_t ro0 = 0.5*(r1+r2);
1487  Double_t fz = (r2-r1)/dz;
1488  Double_t r0sq = point[0]*point[0] + point[1]*point[1];
1489  Double_t rc = ro0 + fz*(point[2]-0.5*(z1+z2));
1490 
1491  Double_t a = dir[0]*dir[0] + dir[1]*dir[1] - fz*fz*dir[2]*dir[2];
1492  Double_t b = point[0]*dir[0] + point[1]*dir[1] - fz*rc*dir[2];
1493  Double_t c = r0sq - rc*rc;
1494 
1495  if (a==0) return TGeoShape::Big();
1496  a = 1./a;
1497  b *= a;
1498  c *= a;
1499  Double_t delta = b*b - c;
1500  if (delta<0) return TGeoShape::Big();
1501  delta = TMath::Sqrt(delta);
1502 
1503  Double_t snxt = -b-delta;
1504  Double_t ptnew[3];
1505  Double_t ddp, phi;
1506  if (snxt>0) {
1507  // check Z range
1508  ptnew[2] = point[2] + snxt*dir[2];
1509  if (((ptnew[2]-z1)*(ptnew[2]-z2)) < 0) {
1510  // check phi range
1511  if (!hasphi) return snxt;
1512  ptnew[0] = point[0] + snxt*dir[0];
1513  ptnew[1] = point[1] + snxt*dir[1];
1514  phi = TMath::ATan2(ptnew[1], ptnew[0]) * TMath::RadToDeg();
1515  if (phi < 0 ) phi+=360.;
1516  ddp = phi-phi1;
1517  if (ddp < 0) ddp+=360.;
1518  // printf("snxt1=%f phi=%f ddp=%f\n", snxt, phi, ddp);
1519  if (ddp<=dphi) return snxt;
1520  }
1521  }
1522  snxt = -b+delta;
1523  if (snxt>0) {
1524  // check Z range
1525  ptnew[2] = point[2] + snxt*dir[2];
1526  if (((ptnew[2]-z1)*(ptnew[2]-z2)) < 0) {
1527  // check phi range
1528  if (!hasphi) return snxt;
1529  ptnew[0] = point[0] + snxt*dir[0];
1530  ptnew[1] = point[1] + snxt*dir[1];
1531  phi = TMath::ATan2(ptnew[1], ptnew[0]) * TMath::RadToDeg();
1532  if (phi < 0 ) phi+=360.;
1533  ddp = phi-phi1;
1534  if (ddp < 0) ddp+=360.;
1535  // printf("snxt2=%f phi=%f ddp=%f\n", snxt, phi, ddp);
1536  if (ddp<=dphi) return snxt;
1537  }
1538  }
1539  return TGeoShape::Big();
1540 }
1541 
1542 ////////////////////////////////////////////////////////////////////////////////
1543 /// compute distance from inside point to surface of the tube segment
1544 
1545 Double_t TGeoConeSeg::DistFromInsideS(const Double_t *point, const Double_t *dir, Double_t dz,
1546  Double_t rmin1, Double_t rmax1, Double_t rmin2, Double_t rmax2,
1547  Double_t c1, Double_t s1, Double_t c2, Double_t s2, Double_t cm, Double_t sm, Double_t cdfi)
1548 {
1549  if (dz<=0) return TGeoShape::Big();
1550  // Do Z
1551  Double_t scone = TGeoCone::DistFromInsideS(point,dir,dz,rmin1,rmax1,rmin2,rmax2);
1552  if (scone<=0) return 0.0;
1553  Double_t sfmin = TGeoShape::Big();
1554  Double_t rsq = point[0]*point[0]+point[1]*point[1];
1555  Double_t r = TMath::Sqrt(rsq);
1556  Double_t cpsi=point[0]*cm+point[1]*sm;
1557  if (cpsi>r*cdfi+TGeoShape::Tolerance()) {
1558  sfmin = TGeoShape::DistToPhiMin(point, dir, s1, c1, s2, c2, sm, cm);
1559  return TMath::Min(scone,sfmin);
1560  }
1561  // Point on the phi boundary or outside
1562  // which one: phi1 or phi2
1563  Double_t ddotn, xi, yi;
1564  if (TMath::Abs(point[1]-s1*r) < TMath::Abs(point[1]-s2*r)) {
1565  ddotn = s1*dir[0]-c1*dir[1];
1566  if (ddotn>=0) return 0.0;
1567  ddotn = -s2*dir[0]+c2*dir[1];
1568  if (ddotn<=0) return scone;
1569  sfmin = s2*point[0]-c2*point[1];
1570  if (sfmin<=0) return scone;
1571  sfmin /= ddotn;
1572  if (sfmin >= scone) return scone;
1573  xi = point[0]+sfmin*dir[0];
1574  yi = point[1]+sfmin*dir[1];
1575  if (yi*cm-xi*sm<0) return scone;
1576  return sfmin;
1577  }
1578  ddotn = -s2*dir[0]+c2*dir[1];
1579  if (ddotn>=0) return 0.0;
1580  ddotn = s1*dir[0]-c1*dir[1];
1581  if (ddotn<=0) return scone;
1582  sfmin = -s1*point[0]+c1*point[1];
1583  if (sfmin<=0) return scone;
1584  sfmin /= ddotn;
1585  if (sfmin >= scone) return scone;
1586  xi = point[0]+sfmin*dir[0];
1587  yi = point[1]+sfmin*dir[1];
1588  if (yi*cm-xi*sm>0) return scone;
1589  return sfmin;
1590 }
1591 
1592 ////////////////////////////////////////////////////////////////////////////////
1593 /// compute distance from inside point to surface of the tube segment
1594 
1595 Double_t TGeoConeSeg::DistFromInside(const Double_t *point, const Double_t *dir, Int_t iact, Double_t step, Double_t *safe) const
1596 {
1597  if (iact<3 && safe) {
1598  *safe = TGeoConeSeg::SafetyS(point, kTRUE, fDz,fRmin1,fRmax1,fRmin2,fRmax2,fPhi1,fPhi2);
1599  if (iact==0) return TGeoShape::Big();
1600  if ((iact==1) && (*safe>step)) return TGeoShape::Big();
1601  }
1602  if ((fPhi2-fPhi1)>=360.) return TGeoCone::DistFromInsideS(point,dir,fDz,fRmin1,fRmax1,fRmin2,fRmax2);
1603 
1604  // compute distance to surface
1605  return TGeoConeSeg::DistFromInsideS(point,dir,fDz,fRmin1,fRmax1,fRmin2,fRmax2,fC1,fS1,fC2,fS2,fCm,fSm,fCdfi);
1606 }
1607 
1608 ////////////////////////////////////////////////////////////////////////////////
1609 /// compute distance from outside point to surface of arbitrary tube
1610 
1611 Double_t TGeoConeSeg::DistFromOutsideS(const Double_t *point, const Double_t *dir, Double_t dz,
1612  Double_t rmin1, Double_t rmax1, Double_t rmin2, Double_t rmax2,
1613  Double_t c1, Double_t s1, Double_t c2, Double_t s2, Double_t cm, Double_t sm, Double_t cdfi)
1614 {
1615  if (dz<=0) return TGeoShape::Big();
1616  Double_t r2, cpsi;
1617  // check Z planes
1618  Double_t xi, yi, zi;
1619  Double_t b,delta;
1620  zi = dz - TMath::Abs(point[2]);
1621  Double_t rin,rout;
1622  Double_t s = TGeoShape::Big();
1623  Double_t snxt=TGeoShape::Big();
1624  Bool_t in = kFALSE;
1625  Bool_t inz = (zi<0)?kFALSE:kTRUE;
1626  if (!inz) {
1627  if (point[2]*dir[2]>=0) return TGeoShape::Big();
1628  s = -zi/TMath::Abs(dir[2]);
1629  xi = point[0]+s*dir[0];
1630  yi = point[1]+s*dir[1];
1631  r2=xi*xi+yi*yi;
1632  if (dir[2]>0) {
1633  rin = rmin1;
1634  rout = rmax1;
1635  } else {
1636  rin = rmin2;
1637  rout = rmax2;
1638  }
1639  if ((rin*rin<=r2) && (r2<=rout*rout)) {
1640  cpsi=xi*cm+yi*sm;
1641  if (cpsi>=(cdfi*TMath::Sqrt(r2))) return s;
1642  }
1643  }
1644  Double_t zinv = 1./dz;
1645  Double_t rsq = point[0]*point[0]+point[1]*point[1];
1646  Double_t r = TMath::Sqrt(rsq);
1647  Double_t ro1=0.5*(rmin1+rmin2);
1648  Bool_t hasrmin = (ro1>0)?kTRUE:kFALSE;
1649  Double_t tg1 = 0.0;
1650  Bool_t inrmin = kFALSE;
1651  rin = 0.0;
1652  if (hasrmin) {
1653  tg1=0.5*(rmin2-rmin1)*zinv;
1654  rin = ro1+tg1*point[2];
1655  if (rsq > rin*(rin-TGeoShape::Tolerance())) inrmin=kTRUE;
1656  } else {
1657  inrmin = kTRUE;
1658  }
1659  Double_t ro2=0.5*(rmax1+rmax2);
1660  Double_t tg2=0.5*(rmax2-rmax1)*zinv;
1661  rout = ro2+tg2*point[2];
1662  Bool_t inrmax = kFALSE;
1663  if (r < rout+TGeoShape::Tolerance()) inrmax = kTRUE;
1664  Bool_t inphi = kFALSE;
1665  cpsi=point[0]*cm+point[1]*sm;
1666  if (cpsi>r*cdfi-TGeoShape::Tolerance()) inphi = kTRUE;
1667  in = inz & inrmin & inrmax & inphi;
1668  // If inside, we are most likely on a boundary within machine precision.
1669  if (in) {
1670  Double_t safphi = (cpsi-r*cdfi)*TMath::Sqrt(1.-cdfi*cdfi);
1671  Double_t safrmin = (hasrmin)?TMath::Abs(r-rin):(TGeoShape::Big());
1672  Double_t safrmax = TMath::Abs(r-rout);
1673  // check if on Z boundaries
1674  if (zi<safrmax && zi<safrmin && zi<safphi) {
1675  if (point[2]*dir[2]<0) return 0.0;
1676  return TGeoShape::Big();
1677  }
1678  // check if on Rmax boundary
1679  if (safrmax<safrmin && safrmax<safphi) {
1680  Double_t ddotn = point[0]*dir[0]+point[1]*dir[1]-tg2*dir[2]*r;
1681  if (ddotn<=0) return 0.0;
1682  return TGeoShape::Big();
1683  }
1684  // check if on phi boundary
1685  if (safphi<safrmin) {
1686  // We may cross again a phi of rmin boundary
1687  // check first if we are on phi1 or phi2
1688  Double_t un;
1689  if (point[0]*c1 + point[1]*s1 > point[0]*c2 + point[1]*s2) {
1690  un = dir[0]*s1-dir[1]*c1;
1691  if (un < 0) return 0.0;
1692  if (cdfi>=0) return TGeoShape::Big();
1693  un = -dir[0]*s2+dir[1]*c2;
1694  if (un<0) {
1695  s = -point[0]*s2+point[1]*c2;
1696  if (s>0) {
1697  s /= (-un);
1698  zi = point[2]+s*dir[2];
1699  if (TMath::Abs(zi)<=dz) {
1700  xi = point[0]+s*dir[0];
1701  yi = point[1]+s*dir[1];
1702  if ((yi*cm-xi*sm)>0) {
1703  r2=xi*xi+yi*yi;
1704  rin = ro1+tg1*zi;
1705  rout = ro2+tg2*zi;
1706  if ((rin*rin<=r2) && (rout*rout>=r2)) return s;
1707  }
1708  }
1709  }
1710  }
1711  } else {
1712  un = -dir[0]*s2+dir[1]*c2;
1713  if (un < 0) return 0.0;
1714  if (cdfi>=0) return TGeoShape::Big();
1715  un = dir[0]*s1-dir[1]*c1;
1716  if (un<0) {
1717  s = point[0]*s1-point[1]*c1;
1718  if (s>0) {
1719  s /= (-un);
1720  zi = point[2]+s*dir[2];
1721  if (TMath::Abs(zi)<=dz) {
1722  xi = point[0]+s*dir[0];
1723  yi = point[1]+s*dir[1];
1724  if ((yi*cm-xi*sm)<0) {
1725  r2=xi*xi+yi*yi;
1726  rin = ro1+tg1*zi;
1727  rout = ro2+tg2*zi;
1728  if ((rin*rin<=r2) && (rout*rout>=r2)) return s;
1729  }
1730  }
1731  }
1732  }
1733  }
1734  // We may also cross rmin, second solution coming from outside
1735  Double_t ddotn = point[0]*dir[0]+point[1]*dir[1]-tg1*dir[2]*r;
1736  if (ddotn>=0) return TGeoShape::Big();
1737  if (cdfi>=0) return TGeoShape::Big();
1738  TGeoCone::DistToCone(point, dir, dz, rmin1, rmin2, b, delta);
1739  if (delta<0) return TGeoShape::Big();
1740  snxt = -b-delta;
1741  if (snxt<0) return TGeoShape::Big();
1742  snxt = -b+delta;
1743  zi = point[2]+snxt*dir[2];
1744  if (TMath::Abs(zi)>dz) return TGeoShape::Big();
1745  xi = point[0]+snxt*dir[0];
1746  yi = point[1]+snxt*dir[1];
1747  r2=xi*xi+yi*yi;
1748  cpsi=xi*cm+yi*sm;
1749  if (cpsi>=(cdfi*TMath::Sqrt(r2))) return snxt;
1750  return TGeoShape::Big();
1751  }
1752  // We are on rmin boundary: we may cross again rmin or a phi facette
1753  Double_t ddotn = point[0]*dir[0]+point[1]*dir[1]-tg1*dir[2]*r;
1754  if (ddotn>=0) return 0.0;
1755  TGeoCone::DistToCone(point, dir, dz, rmin1, rmin2, b, delta);
1756  if (delta<0) return 0.0;
1757  snxt = -b+delta;
1758  if (snxt<0) return TGeoShape::Big();
1759  if (TMath::Abs(-b-delta)>snxt) return TGeoShape::Big();
1760  zi = point[2]+snxt*dir[2];
1761  if (TMath::Abs(zi)>dz) return TGeoShape::Big();
1762  // OK, we cross rmin at snxt - check if within phi range
1763  xi = point[0]+snxt*dir[0];
1764  yi = point[1]+snxt*dir[1];
1765  r2=xi*xi+yi*yi;
1766  cpsi=xi*cm+yi*sm;
1767  if (cpsi>=(cdfi*TMath::Sqrt(r2))) return snxt;
1768  // we cross rmin in the phi gap - we may cross a phi facette
1769  if (cdfi>=0) return TGeoShape::Big();
1770  Double_t un=-dir[0]*s1+dir[1]*c1;
1771  if (un > 0) {
1772  s=point[0]*s1-point[1]*c1;
1773  if (s>=0) {
1774  s /= un;
1775  zi=point[2]+s*dir[2];
1776  if (TMath::Abs(zi)<=dz) {
1777  xi=point[0]+s*dir[0];
1778  yi=point[1]+s*dir[1];
1779  if ((yi*cm-xi*sm)<=0) {
1780  r2=xi*xi+yi*yi;
1781  rin = ro1+tg1*zi;
1782  rout = ro2+tg2*zi;
1783  if ((rin*rin<=r2) && (rout*rout>=r2)) return s;
1784  }
1785  }
1786  }
1787  }
1788  un=dir[0]*s2-dir[1]*c2;
1789  if (un > 0) {
1790  s=(point[1]*c2-point[0]*s2)/un;
1791  if (s>=0) {
1792  zi=point[2]+s*dir[2];
1793  if (TMath::Abs(zi)<=dz) {
1794  xi=point[0]+s*dir[0];
1795  yi=point[1]+s*dir[1];
1796  if ((yi*cm-xi*sm)>=0) {
1797  r2=xi*xi+yi*yi;
1798  rin = ro1+tg1*zi;
1799  rout = ro2+tg2*zi;
1800  if ((rin*rin<=r2) && (rout*rout>=r2)) return s;
1801  }
1802  }
1803  }
1804  }
1805  return TGeoShape::Big();
1806  }
1807 
1808  // The point is really outside
1809  Double_t sr1 = TGeoShape::Big();
1810  if (!inrmax) {
1811  // check crossing with outer cone
1812  TGeoCone::DistToCone(point, dir, dz, rmax1, rmax2, b, delta);
1813  if (delta>=0) {
1814  s = -b-delta;
1815  if (s>0) {
1816  zi=point[2]+s*dir[2];
1817  if (TMath::Abs(zi)<=dz) {
1818  xi=point[0]+s*dir[0];
1819  yi=point[1]+s*dir[1];
1820  r2=xi*xi+yi*yi;
1821  cpsi=xi*cm+yi*sm;
1822  if (cpsi>=(cdfi*TMath::Sqrt(r2))) return s; // rmax crossing
1823  }
1824  }
1825  s = -b+delta;
1826  if (s>0) {
1827  zi=point[2]+s*dir[2];
1828  if (TMath::Abs(zi)<=dz) {
1829  xi=point[0]+s*dir[0];
1830  yi=point[1]+s*dir[1];
1831  r2=xi*xi+yi*yi;
1832  cpsi=xi*cm+yi*sm;
1833  if (cpsi>=(cdfi*TMath::Sqrt(r2))) sr1=s;
1834  }
1835  }
1836  }
1837  }
1838  // check crossing with inner cone
1839  Double_t sr2 = TGeoShape::Big();
1840  TGeoCone::DistToCone(point, dir, dz, rmin1, rmin2, b, delta);
1841  if (delta>=0) {
1842  s = -b-delta;
1843  if (s>0) {
1844  zi=point[2]+s*dir[2];
1845  if (TMath::Abs(zi)<=dz) {
1846  xi=point[0]+s*dir[0];
1847  yi=point[1]+s*dir[1];
1848  r2=xi*xi+yi*yi;
1849  cpsi=xi*cm+yi*sm;
1850  if (cpsi>=(cdfi*TMath::Sqrt(r2))) sr2=s;
1851  }
1852  }
1853  if (sr2>1E10) {
1854  s = -b+delta;
1855  if (s>0) {
1856  zi=point[2]+s*dir[2];
1857  if (TMath::Abs(zi)<=dz) {
1858  xi=point[0]+s*dir[0];
1859  yi=point[1]+s*dir[1];
1860  r2=xi*xi+yi*yi;
1861  cpsi=xi*cm+yi*sm;
1862  if (cpsi>=(cdfi*TMath::Sqrt(r2))) sr2=s;
1863  }
1864  }
1865  }
1866  }
1867  snxt = TMath::Min(sr1,sr2);
1868  // Check phi crossing
1869  s = TGeoShape::DistToPhiMin(point,dir,s1,c1,s2,c2,sm,cm,kFALSE);
1870  if (s>snxt) return snxt;
1871  zi=point[2]+s*dir[2];
1872  if (TMath::Abs(zi)>dz) return snxt;
1873  xi=point[0]+s*dir[0];
1874  yi=point[1]+s*dir[1];
1875  r2=xi*xi+yi*yi;
1876  rout = ro2+tg2*zi;
1877  if (r2>rout*rout) return snxt;
1878  rin = ro1+tg1*zi;
1879  if (r2>=rin*rin) return s; // phi crossing
1880  return snxt;
1881 }
1882 
1883 ////////////////////////////////////////////////////////////////////////////////
1884 /// compute distance from outside point to surface of the tube
1885 
1886 Double_t TGeoConeSeg::DistFromOutside(const Double_t *point, const Double_t *dir, Int_t iact, Double_t step, Double_t *safe) const
1887 {
1888  // compute safe radius
1889  if (iact<3 && safe) {
1890  *safe = Safety(point, kFALSE);
1891  if (iact==0) return TGeoShape::Big();
1892  if ((iact==1) && (*safe>step)) return TGeoShape::Big();
1893  }
1894  // Check if the bounding box is crossed within the requested distance
1895  Double_t sdist = TGeoBBox::DistFromOutside(point,dir, fDX, fDY, fDZ, fOrigin, step);
1896  if (sdist>=step) return TGeoShape::Big();
1897  if ((fPhi2-fPhi1)>=360.) return TGeoCone::DistFromOutsideS(point,dir,fDz,fRmin1,fRmax1,fRmin2,fRmax2);
1898  return TGeoConeSeg::DistFromOutsideS(point,dir,fDz,fRmin1,fRmax1,fRmin2,fRmax2,fC1,fS1,fC2,fS2,fCm,fSm,fCdfi);
1899 }
1900 
1901 ////////////////////////////////////////////////////////////////////////////////
1902 /// compute closest distance from point px,py to each corner
1903 
1904 Int_t TGeoConeSeg::DistancetoPrimitive(Int_t px, Int_t py)
1905 {
1906  Int_t n = gGeoManager->GetNsegments()+1;
1907  const Int_t numPoints = 4*n;
1908  return ShapeDistancetoPrimitive(numPoints, px, py);
1909 }
1910 
1911 ////////////////////////////////////////////////////////////////////////////////
1912 /// Divide this cone segment shape belonging to volume "voldiv" into ndiv volumes
1913 /// called divname, from start position with the given step. Returns pointer
1914 /// to created division cell volume in case of Z divisions. For Z division
1915 /// creates all volumes with different shapes and returns pointer to volume that
1916 /// was divided. In case a wrong division axis is supplied, returns pointer to
1917 /// volume that was divided.
1918 
1919 TGeoVolume *TGeoConeSeg::Divide(TGeoVolume *voldiv, const char *divname, Int_t iaxis, Int_t ndiv,
1920  Double_t start, Double_t step)
1921 {
1922  TGeoShape *shape; //--- shape to be created
1923  TGeoVolume *vol; //--- division volume to be created
1924  TGeoVolumeMulti *vmulti; //--- generic divided volume
1925  TGeoPatternFinder *finder; //--- finder to be attached
1926  TString opt = ""; //--- option to be attached
1927  Double_t dphi;
1928  Int_t id;
1929  Double_t end = start+ndiv*step;
1930  switch (iaxis) {
1931  case 1: //--- R division
1932  Error("Divide","division of a cone segment on R not implemented");
1933  return 0;
1934  case 2: //--- Phi division
1935  dphi = fPhi2-fPhi1;
1936  if (dphi<0) dphi+=360.;
1937  finder = new TGeoPatternCylPhi(voldiv, ndiv, start, end);
1938  voldiv->SetFinder(finder);
1939  finder->SetDivIndex(voldiv->GetNdaughters());
1940  shape = new TGeoConeSeg(fDz, fRmin1, fRmax1, fRmin2, fRmax2, -step/2, step/2);
1941  vol = new TGeoVolume(divname, shape, voldiv->GetMedium());
1942  vmulti = gGeoManager->MakeVolumeMulti(divname, voldiv->GetMedium());
1943  vmulti->AddVolume(vol);
1944  opt = "Phi";
1945  for (id=0; id<ndiv; id++) {
1946  voldiv->AddNodeOffset(vol, id, start+id*step+step/2, opt.Data());
1947  ((TGeoNodeOffset*)voldiv->GetNodes()->At(voldiv->GetNdaughters()-1))->SetFinder(finder);
1948  }
1949  return vmulti;
1950  case 3: //--- Z division
1951  finder = new TGeoPatternZ(voldiv, ndiv, start, end);
1952  vmulti = gGeoManager->MakeVolumeMulti(divname, voldiv->GetMedium());
1953  voldiv->SetFinder(finder);
1954  finder->SetDivIndex(voldiv->GetNdaughters());
1955  for (id=0; id<ndiv; id++) {
1956  Double_t z1 = start+id*step;
1957  Double_t z2 = start+(id+1)*step;
1958  Double_t rmin1n = 0.5*(fRmin1*(fDz-z1)+fRmin2*(fDz+z1))/fDz;
1959  Double_t rmax1n = 0.5*(fRmax1*(fDz-z1)+fRmax2*(fDz+z1))/fDz;
1960  Double_t rmin2n = 0.5*(fRmin1*(fDz-z2)+fRmin2*(fDz+z2))/fDz;
1961  Double_t rmax2n = 0.5*(fRmax1*(fDz-z2)+fRmax2*(fDz+z2))/fDz;
1962  shape = new TGeoConeSeg(step/2, rmin1n, rmax1n, rmin2n, rmax2n, fPhi1, fPhi2);
1963  vol = new TGeoVolume(divname, shape, voldiv->GetMedium());
1964  vmulti->AddVolume(vol);
1965  opt = "Z";
1966  voldiv->AddNodeOffset(vol, id, start+id*step+step/2, opt.Data());
1967  ((TGeoNodeOffset*)voldiv->GetNodes()->At(voldiv->GetNdaughters()-1))->SetFinder(finder);
1968  }
1969  return vmulti;
1970  default:
1971  Error("Divide", "Wrong axis type for division");
1972  return 0;
1973  }
1974 }
1975 
1976 ////////////////////////////////////////////////////////////////////////////////
1977 /// Get range of shape for a given axis.
1978 
1979 Double_t TGeoConeSeg::GetAxisRange(Int_t iaxis, Double_t &xlo, Double_t &xhi) const
1980 {
1981  xlo = 0;
1982  xhi = 0;
1983  Double_t dx = 0;
1984  switch (iaxis) {
1985  case 2:
1986  xlo = fPhi1;
1987  xhi = fPhi2;
1988  dx = xhi-xlo;
1989  return dx;
1990  case 3:
1991  xlo = -fDz;
1992  xhi = fDz;
1993  dx = xhi-xlo;
1994  return dx;
1995  }
1996  return dx;
1997 }
1998 
1999 ////////////////////////////////////////////////////////////////////////////////
2000 /// Fill vector param[4] with the bounding cylinder parameters. The order
2001 /// is the following : Rmin, Rmax, Phi1, Phi2
2002 
2003 void TGeoConeSeg::GetBoundingCylinder(Double_t *param) const
2004 {
2005  param[0] = TMath::Min(fRmin1, fRmin2); // Rmin
2006  param[0] *= param[0];
2007  param[1] = TMath::Max(fRmax1, fRmax2); // Rmax
2008  param[1] *= param[1];
2009  param[2] = (fPhi1<0)?(fPhi1+360.):fPhi1; // Phi1
2010  param[3] = fPhi2; // Phi2
2011  while (param[3]<param[2]) param[3]+=360.;
2012 }
2013 
2014 ////////////////////////////////////////////////////////////////////////////////
2015 /// in case shape has some negative parameters, these has to be computed
2016 /// in order to fit the mother
2017 
2018 TGeoShape *TGeoConeSeg::GetMakeRuntimeShape(TGeoShape *mother, TGeoMatrix * /*mat*/) const
2019 {
2020  if (!TestShapeBit(kGeoRunTimeShape)) return 0;
2021  if (!mother->TestShapeBit(kGeoConeSeg)) {
2022  Error("GetMakeRuntimeShape", "invalid mother");
2023  return 0;
2024  }
2025  Double_t rmin1, rmax1, rmin2, rmax2, dz;
2026  rmin1 = fRmin1;
2027  rmax1 = fRmax1;
2028  rmin2 = fRmin2;
2029  rmax2 = fRmax2;
2030  dz = fDz;
2031  if (fDz<0) dz=((TGeoCone*)mother)->GetDz();
2032  if (fRmin1<0)
2033  rmin1 = ((TGeoCone*)mother)->GetRmin1();
2034  if ((fRmax1<0) || (fRmax1<fRmin1))
2035  rmax1 = ((TGeoCone*)mother)->GetRmax1();
2036  if (fRmin2<0)
2037  rmin2 = ((TGeoCone*)mother)->GetRmin2();
2038  if ((fRmax2<0) || (fRmax2<fRmin2))
2039  rmax2 = ((TGeoCone*)mother)->GetRmax2();
2040 
2041  return (new TGeoConeSeg(GetName(), dz, rmin1, rmax1, rmin2, rmax2, fPhi1, fPhi2));
2042 }
2043 
2044 ////////////////////////////////////////////////////////////////////////////////
2045 /// print shape parameters
2046 
2047 void TGeoConeSeg::InspectShape() const
2048 {
2049  printf("*** Shape %s: TGeoConeSeg ***\n", GetName());
2050  printf(" dz = %11.5f\n", fDz);
2051  printf(" Rmin1 = %11.5f\n", fRmin1);
2052  printf(" Rmax1 = %11.5f\n", fRmax1);
2053  printf(" Rmin2 = %11.5f\n", fRmin2);
2054  printf(" Rmax2 = %11.5f\n", fRmax2);
2055  printf(" phi1 = %11.5f\n", fPhi1);
2056  printf(" phi2 = %11.5f\n", fPhi2);
2057  printf(" Bounding box:\n");
2058  TGeoBBox::InspectShape();
2059 }
2060 
2061  ///////////////////////////////////////////////////////////////////////////////
2062  /// Creates a TBuffer3D describing *this* shape.
2063  /// Coordinates are in local reference frame.
2064 
2065 TBuffer3D *TGeoConeSeg::MakeBuffer3D() const
2066 {
2067  Int_t n = gGeoManager->GetNsegments()+1;
2068  Int_t nbPnts = 4*n;
2069  Int_t nbSegs = 2*nbPnts;
2070  Int_t nbPols = nbPnts-2;
2071  TBuffer3D* buff = new TBuffer3D(TBuffer3DTypes::kGeneric,
2072  nbPnts, 3*nbPnts,
2073  nbSegs, 3*nbSegs,
2074  nbPols, 6*nbPols);
2075 
2076  if (buff)
2077  {
2078  SetPoints(buff->fPnts);
2079  SetSegsAndPols(*buff);
2080  }
2081 
2082  return buff;
2083 }
2084 
2085 ////////////////////////////////////////////////////////////////////////////////
2086 /// Fill TBuffer3D structure for segments and polygons.
2087 
2088 void TGeoConeSeg::SetSegsAndPols(TBuffer3D &buffer) const
2089 {
2090  Int_t i, j;
2091  Int_t n = gGeoManager->GetNsegments()+1;
2092  Int_t c = GetBasicColor();
2093 
2094  memset(buffer.fSegs, 0, buffer.NbSegs()*3*sizeof(Int_t));
2095  for (i = 0; i < 4; i++) {
2096  for (j = 1; j < n; j++) {
2097  buffer.fSegs[(i*n+j-1)*3 ] = c;
2098  buffer.fSegs[(i*n+j-1)*3+1] = i*n+j-1;
2099  buffer.fSegs[(i*n+j-1)*3+2] = i*n+j;
2100  }
2101  }
2102  for (i = 4; i < 6; i++) {
2103  for (j = 0; j < n; j++) {
2104  buffer.fSegs[(i*n+j)*3 ] = c+1;
2105  buffer.fSegs[(i*n+j)*3+1] = (i-4)*n+j;
2106  buffer.fSegs[(i*n+j)*3+2] = (i-2)*n+j;
2107  }
2108  }
2109  for (i = 6; i < 8; i++) {
2110  for (j = 0; j < n; j++) {
2111  buffer.fSegs[(i*n+j)*3 ] = c;
2112  buffer.fSegs[(i*n+j)*3+1] = 2*(i-6)*n+j;
2113  buffer.fSegs[(i*n+j)*3+2] = (2*(i-6)+1)*n+j;
2114  }
2115  }
2116 
2117  Int_t indx = 0;
2118  memset(buffer.fPols, 0, buffer.NbPols()*6*sizeof(Int_t));
2119  i = 0;
2120  for (j = 0; j < n-1; j++) {
2121  buffer.fPols[indx++] = c;
2122  buffer.fPols[indx++] = 4;
2123  buffer.fPols[indx++] = (4+i)*n+j+1;
2124  buffer.fPols[indx++] = (2+i)*n+j;
2125  buffer.fPols[indx++] = (4+i)*n+j;
2126  buffer.fPols[indx++] = i*n+j;
2127  }
2128  i = 1;
2129  for (j = 0; j < n-1; j++) {
2130  buffer.fPols[indx++] = c;
2131  buffer.fPols[indx++] = 4;
2132  buffer.fPols[indx++] = i*n+j;
2133  buffer.fPols[indx++] = (4+i)*n+j;
2134  buffer.fPols[indx++] = (2+i)*n+j;
2135  buffer.fPols[indx++] = (4+i)*n+j+1;
2136  }
2137  i = 2;
2138  for (j = 0; j < n-1; j++) {
2139  buffer.fPols[indx++] = c+i;
2140  buffer.fPols[indx++] = 4;
2141  buffer.fPols[indx++] = (i-2)*2*n+j;
2142  buffer.fPols[indx++] = (4+i)*n+j;
2143  buffer.fPols[indx++] = ((i-2)*2+1)*n+j;
2144  buffer.fPols[indx++] = (4+i)*n+j+1;
2145  }
2146  i = 3;
2147  for (j = 0; j < n-1; j++) {
2148  buffer.fPols[indx++] = c+i;
2149  buffer.fPols[indx++] = 4;
2150  buffer.fPols[indx++] = (4+i)*n+j+1;
2151  buffer.fPols[indx++] = ((i-2)*2+1)*n+j;
2152  buffer.fPols[indx++] = (4+i)*n+j;
2153  buffer.fPols[indx++] = (i-2)*2*n+j;
2154  }
2155  buffer.fPols[indx++] = c+2;
2156  buffer.fPols[indx++] = 4;
2157  buffer.fPols[indx++] = 6*n;
2158  buffer.fPols[indx++] = 4*n;
2159  buffer.fPols[indx++] = 7*n;
2160  buffer.fPols[indx++] = 5*n;
2161  buffer.fPols[indx++] = c+2;
2162  buffer.fPols[indx++] = 4;
2163  buffer.fPols[indx++] = 6*n-1;
2164  buffer.fPols[indx++] = 8*n-1;
2165  buffer.fPols[indx++] = 5*n-1;
2166  buffer.fPols[indx++] = 7*n-1;
2167 }
2168 
2169 ////////////////////////////////////////////////////////////////////////////////
2170 /// computes the closest distance from given point to this shape, according
2171 /// to option. The matching point on the shape is stored in spoint.
2172 
2173 Double_t TGeoConeSeg::Safety(const Double_t *point, Bool_t in) const
2174 {
2175  Double_t safe = TGeoCone::Safety(point,in);
2176  if ((fPhi2-fPhi1)>=360.) return safe;
2177  Double_t safphi = TGeoShape::SafetyPhi(point, in, fPhi1, fPhi2);
2178  if (in) return TMath::Min(safe, safphi);
2179  if (safe>1.E10) return safphi;
2180  return TMath::Max(safe, safphi);
2181 }
2182 
2183 ////////////////////////////////////////////////////////////////////////////////
2184 /// Static method to compute the closest distance from given point to this shape.
2185 
2186 Double_t TGeoConeSeg::SafetyS(const Double_t *point, Bool_t in, Double_t dz, Double_t rmin1, Double_t rmax1,
2187  Double_t rmin2, Double_t rmax2, Double_t phi1, Double_t phi2, Int_t skipz)
2188 {
2189  Double_t safe = TGeoCone::SafetyS(point,in,dz,rmin1,rmax1,rmin2,rmax2,skipz);
2190  if ((phi2-phi1)>=360.) return safe;
2191  Double_t safphi = TGeoShape::SafetyPhi(point,in,phi1,phi2);
2192  if (in) return TMath::Min(safe, safphi);
2193  if (safe>1.E10) return safphi;
2194  return TMath::Max(safe, safphi);
2195 }
2196 
2197 ////////////////////////////////////////////////////////////////////////////////
2198 /// Save a primitive as a C++ statement(s) on output stream "out".
2199 
2200 void TGeoConeSeg::SavePrimitive(std::ostream &out, Option_t * /*option*/ /*= ""*/)
2201 {
2202  if (TObject::TestBit(kGeoSavePrimitive)) return;
2203  out << " // Shape: " << GetName() << " type: " << ClassName() << std::endl;
2204  out << " dz = " << fDz << ";" << std::endl;
2205  out << " rmin1 = " << fRmin1 << ";" << std::endl;
2206  out << " rmax1 = " << fRmax1 << ";" << std::endl;
2207  out << " rmin2 = " << fRmin2 << ";" << std::endl;
2208  out << " rmax2 = " << fRmax2 << ";" << std::endl;
2209  out << " phi1 = " << fPhi1 << ";" << std::endl;
2210  out << " phi2 = " << fPhi2 << ";" << std::endl;
2211  out << " TGeoShape *" << GetPointerName() << " = new TGeoConeSeg(\"" << GetName() << "\", dz,rmin1,rmax1,rmin2,rmax2,phi1,phi2);" << std::endl;
2212  TObject::SetBit(TGeoShape::kGeoSavePrimitive);
2213 }
2214 
2215 ////////////////////////////////////////////////////////////////////////////////
2216 /// Set dimensions of the cone segment.
2217 
2218 void TGeoConeSeg::SetConsDimensions(Double_t dz, Double_t rmin1, Double_t rmax1,
2219  Double_t rmin2, Double_t rmax2, Double_t phi1, Double_t phi2)
2220 {
2221  fDz = dz;
2222  fRmin1 = rmin1;
2223  fRmax1 = rmax1;
2224  fRmin2 = rmin2;
2225  fRmax2 = rmax2;
2226  fPhi1 = phi1;
2227  while (fPhi1<0) fPhi1+=360.;
2228  fPhi2 = phi2;
2229  while (fPhi2<=fPhi1) fPhi2+=360.;
2230  if (TGeoShape::IsSameWithinTolerance(fPhi1,fPhi2)) Fatal("SetConsDimensions", "In shape %s invalid phi1=%g, phi2=%g\n", GetName(), fPhi1, fPhi2);
2231  InitTrigonometry();
2232 }
2233 
2234 ////////////////////////////////////////////////////////////////////////////////
2235 /// Set dimensions of the cone segment from an array.
2236 
2237 void TGeoConeSeg::SetDimensions(Double_t *param)
2238 {
2239  Double_t dz = param[0];
2240  Double_t rmin1 = param[1];
2241  Double_t rmax1 = param[2];
2242  Double_t rmin2 = param[3];
2243  Double_t rmax2 = param[4];
2244  Double_t phi1 = param[5];
2245  Double_t phi2 = param[6];
2246  SetConsDimensions(dz, rmin1, rmax1,rmin2, rmax2, phi1, phi2);
2247 }
2248 
2249 ////////////////////////////////////////////////////////////////////////////////
2250 /// Create cone segment mesh points.
2251 
2252 void TGeoConeSeg::SetPoints(Double_t *points) const
2253 {
2254  Int_t j, n;
2255  Float_t dphi,phi,phi1, phi2,dz;
2256 
2257  n = gGeoManager->GetNsegments()+1;
2258  dz = fDz;
2259  phi1 = fPhi1;
2260  phi2 = fPhi2;
2261 
2262  dphi = (phi2-phi1)/(n-1);
2263 
2264  Int_t indx = 0;
2265 
2266  if (points) {
2267  for (j = 0; j < n; j++) {
2268  phi = (fPhi1+j*dphi)*TMath::DegToRad();
2269  points[indx++] = fRmin1 * TMath::Cos(phi);
2270  points[indx++] = fRmin1 * TMath::Sin(phi);
2271  points[indx++] = -dz;
2272  }
2273  for (j = 0; j < n; j++) {
2274  phi = (fPhi1+j*dphi)*TMath::DegToRad();
2275  points[indx++] = fRmax1 * TMath::Cos(phi);
2276  points[indx++] = fRmax1 * TMath::Sin(phi);
2277  points[indx++] = -dz;
2278  }
2279  for (j = 0; j < n; j++) {
2280  phi = (fPhi1+j*dphi)*TMath::DegToRad();
2281  points[indx++] = fRmin2 * TMath::Cos(phi);
2282  points[indx++] = fRmin2 * TMath::Sin(phi);
2283  points[indx++] = dz;
2284  }
2285  for (j = 0; j < n; j++) {
2286  phi = (fPhi1+j*dphi)*TMath::DegToRad();
2287  points[indx++] = fRmax2 * TMath::Cos(phi);
2288  points[indx++] = fRmax2 * TMath::Sin(phi);
2289  points[indx++] = dz;
2290  }
2291  }
2292 }
2293 
2294 ////////////////////////////////////////////////////////////////////////////////
2295 /// Create cone segment mesh points.
2296 
2297 void TGeoConeSeg::SetPoints(Float_t *points) const
2298 {
2299  Int_t j, n;
2300  Float_t dphi,phi,phi1, phi2,dz;
2301 
2302  n = gGeoManager->GetNsegments()+1;
2303  dz = fDz;
2304  phi1 = fPhi1;
2305  phi2 = fPhi2;
2306 
2307  dphi = (phi2-phi1)/(n-1);
2308 
2309  Int_t indx = 0;
2310 
2311  if (points) {
2312  for (j = 0; j < n; j++) {
2313  phi = (fPhi1+j*dphi)*TMath::DegToRad();
2314  points[indx++] = fRmin1 * TMath::Cos(phi);
2315  points[indx++] = fRmin1 * TMath::Sin(phi);
2316  points[indx++] = -dz;
2317  }
2318  for (j = 0; j < n; j++) {
2319  phi = (fPhi1+j*dphi)*TMath::DegToRad();
2320  points[indx++] = fRmax1 * TMath::Cos(phi);
2321  points[indx++] = fRmax1 * TMath::Sin(phi);
2322  points[indx++] = -dz;
2323  }
2324  for (j = 0; j < n; j++) {
2325  phi = (fPhi1+j*dphi)*TMath::DegToRad();
2326  points[indx++] = fRmin2 * TMath::Cos(phi);
2327  points[indx++] = fRmin2 * TMath::Sin(phi);
2328  points[indx++] = dz;
2329  }
2330  for (j = 0; j < n; j++) {
2331  phi = (fPhi1+j*dphi)*TMath::DegToRad();
2332  points[indx++] = fRmax2 * TMath::Cos(phi);
2333  points[indx++] = fRmax2 * TMath::Sin(phi);
2334  points[indx++] = dz;
2335  }
2336  }
2337 }
2338 
2339 ////////////////////////////////////////////////////////////////////////////////
2340 /// Returns numbers of vertices, segments and polygons composing the shape mesh.
2341 
2342 void TGeoConeSeg::GetMeshNumbers(Int_t &nvert, Int_t &nsegs, Int_t &npols) const
2343 {
2344  Int_t n = gGeoManager->GetNsegments()+1;
2345  nvert = n*4;
2346  nsegs = n*8;
2347  npols = n*4-2;
2348 }
2349 
2350 ////////////////////////////////////////////////////////////////////////////////
2351 /// Return number of vertices of the mesh representation
2352 
2353 Int_t TGeoConeSeg::GetNmeshVertices() const
2354 {
2355  Int_t n = gGeoManager->GetNsegments()+1;
2356  Int_t numPoints = n*4;
2357  return numPoints;
2358 }
2359 
2360 ////////////////////////////////////////////////////////////////////////////////
2361 /// Fill size of this 3-D object
2362 
2363 void TGeoConeSeg::Sizeof3D() const
2364 {
2365 }
2366 
2367 ////////////////////////////////////////////////////////////////////////////////
2368 /// Fills a static 3D buffer and returns a reference.
2369 
2370 const TBuffer3D & TGeoConeSeg::GetBuffer3D(Int_t reqSections, Bool_t localFrame) const
2371 {
2372  static TBuffer3D buffer(TBuffer3DTypes::kGeneric);
2373 
2374  TGeoBBox::FillBuffer3D(buffer, reqSections, localFrame);
2375 
2376  if (reqSections & TBuffer3D::kRawSizes) {
2377  Int_t n = gGeoManager->GetNsegments()+1;
2378  Int_t nbPnts = 4*n;
2379  Int_t nbSegs = 2*nbPnts;
2380  Int_t nbPols = nbPnts-2;
2381  if (buffer.SetRawSizes(nbPnts, 3*nbPnts, nbSegs, 3*nbSegs, nbPols, 6*nbPols)) {
2382  buffer.SetSectionsValid(TBuffer3D::kRawSizes);
2383  }
2384  }
2385  if ((reqSections & TBuffer3D::kRaw) && buffer.SectionsValid(TBuffer3D::kRawSizes)) {
2386  SetPoints(buffer.fPnts);
2387  if (!buffer.fLocalFrame) {
2388  TransformPoints(buffer.fPnts, buffer.NbPnts());
2389  }
2390 
2391  SetSegsAndPols(buffer);
2392  buffer.SetSectionsValid(TBuffer3D::kRaw);
2393  }
2394 
2395  return buffer;
2396 }
2397 
2398 ////////////////////////////////////////////////////////////////////////////////
2399 /// Fills array with n random points located on the line segments of the shape mesh.
2400 /// The output array must be provided with a length of minimum 3*npoints. Returns
2401 /// true if operation is implemented.
2402 
2403 Bool_t TGeoConeSeg::GetPointsOnSegments(Int_t npoints, Double_t *array) const
2404 {
2405  if (npoints > (npoints/2)*2) {
2406  Error("GetPointsOnSegments","Npoints must be even number");
2407  return kFALSE;
2408  }
2409  Int_t nc = (Int_t)TMath::Sqrt(0.5*npoints);
2410  Double_t dphi = (fPhi2-fPhi1)*TMath::DegToRad()/(nc-1);
2411  Double_t phi = 0;
2412  Double_t phi1 = fPhi1 * TMath::DegToRad();
2413  Int_t ntop = npoints/2 - nc*(nc-1);
2414  Double_t dz = 2*fDz/(nc-1);
2415  Double_t z = 0;
2416  Double_t rmin = 0.;
2417  Double_t rmax = 0.;
2418  Int_t icrt = 0;
2419  Int_t nphi = nc;
2420  // loop z sections
2421  for (Int_t i=0; i<nc; i++) {
2422  if (i == (nc-1)) {
2423  nphi = ntop;
2424  dphi = (fPhi2-fPhi1)*TMath::DegToRad()/(nphi-1);
2425  }
2426  z = -fDz + i*dz;
2427  rmin = 0.5*(fRmin1+fRmin2) + 0.5*(fRmin2-fRmin1)*z/fDz;
2428  rmax = 0.5*(fRmax1+fRmax2) + 0.5*(fRmax2-fRmax1)*z/fDz;
2429  // loop points on circle sections
2430  for (Int_t j=0; j<nphi; j++) {
2431  phi = phi1 + j*dphi;
2432  array[icrt++] = rmin * TMath::Cos(phi);
2433  array[icrt++] = rmin * TMath::Sin(phi);
2434  array[icrt++] = z;
2435  array[icrt++] = rmax * TMath::Cos(phi);
2436  array[icrt++] = rmax * TMath::Sin(phi);
2437  array[icrt++] = z;
2438  }
2439  }
2440  return kTRUE;
2441 }
2442 
2443 ////////////////////////////////////////////////////////////////////////////////
2444 /// Check the inside status for each of the points in the array.
2445 /// Input: Array of point coordinates + vector size
2446 /// Output: Array of Booleans for the inside of each point
2447 
2448 void TGeoConeSeg::Contains_v(const Double_t *points, Bool_t *inside, Int_t vecsize) const
2449 {
2450  for (Int_t i=0; i<vecsize; i++) inside[i] = Contains(&points[3*i]);
2451 }
2452 
2453 ////////////////////////////////////////////////////////////////////////////////
2454 /// Compute the normal for an array o points so that norm.dot.dir is positive
2455 /// Input: Arrays of point coordinates and directions + vector size
2456 /// Output: Array of normal directions
2457 
2458 void TGeoConeSeg::ComputeNormal_v(const Double_t *points, const Double_t *dirs, Double_t *norms, Int_t vecsize)
2459 {
2460  for (Int_t i=0; i<vecsize; i++) ComputeNormal(&points[3*i], &dirs[3*i], &norms[3*i]);
2461 }
2462 
2463 ////////////////////////////////////////////////////////////////////////////////
2464 /// Compute distance from array of input points having directions specified by dirs. Store output in dists
2465 
2466 void TGeoConeSeg::DistFromInside_v(const Double_t *points, const Double_t *dirs, Double_t *dists, Int_t vecsize, Double_t* step) const
2467 {
2468  for (Int_t i=0; i<vecsize; i++) dists[i] = DistFromInside(&points[3*i], &dirs[3*i], 3, step[i]);
2469 }
2470 
2471 ////////////////////////////////////////////////////////////////////////////////
2472 /// Compute distance from array of input points having directions specified by dirs. Store output in dists
2473 
2474 void TGeoConeSeg::DistFromOutside_v(const Double_t *points, const Double_t *dirs, Double_t *dists, Int_t vecsize, Double_t* step) const
2475 {
2476  for (Int_t i=0; i<vecsize; i++) dists[i] = DistFromOutside(&points[3*i], &dirs[3*i], 3, step[i]);
2477 }
2478 
2479 ////////////////////////////////////////////////////////////////////////////////
2480 /// Compute safe distance from each of the points in the input array.
2481 /// Input: Array of point coordinates, array of statuses for these points, size of the arrays
2482 /// Output: Safety values
2483 
2484 void TGeoConeSeg::Safety_v(const Double_t *points, const Bool_t *inside, Double_t *safe, Int_t vecsize) const
2485 {
2486  for (Int_t i=0; i<vecsize; i++) safe[i] = Safety(&points[3*i], inside[i]);
2487 }