Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
TGeoPcon.cxx
Go to the documentation of this file.
1 // @(#)root/geom:$Id$
2 // Author: Andrei Gheata 24/10/01
3 // TGeoPcon::Contains() 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 TGeoPcon
14 \ingroup Geometry_classes
15 A polycone.
16 
17 It has at least 9 parameters:
18  - the lower phi limit;
19  - the range in phi;
20  - the number of z planes (at least two) where the inner/outer
21  radii are changing;
22  - z coordinate, inner and outer radius for each z plane
23 
24 Begin_Macro(source)
25 {
26  TCanvas *c = new TCanvas("c", "c",0,0,600,600);
27  new TGeoManager("pcon", "poza10");
28  TGeoMaterial *mat = new TGeoMaterial("Al", 26.98,13,2.7);
29  TGeoMedium *med = new TGeoMedium("MED",1,mat);
30  TGeoVolume *top = gGeoManager->MakeBox("TOP",med,100,100,100);
31  gGeoManager->SetTopVolume(top);
32  TGeoVolume *vol = gGeoManager->MakePcon("PCON",med, -30.0,300,4);
33  TGeoPcon *pcon = (TGeoPcon*)(vol->GetShape());
34  pcon->DefineSection(0,0,15,20);
35  pcon->DefineSection(1,20,15,20);
36  pcon->DefineSection(2,20,15,25);
37  pcon->DefineSection(3,50,15,20);
38  vol->SetLineWidth(2);
39  top->AddNode(vol,1);
40  gGeoManager->CloseGeometry();
41  gGeoManager->SetNsegments(30);
42  top->Draw();
43  TView *view = gPad->GetView();
44  view->ShowAxis();
45 }
46 End_Macro
47 */
48 
49 #include "TGeoPcon.h"
50 
51 #include "Riostream.h"
52 #include "TBuffer.h"
53 #include "TGeoManager.h"
54 #include "TGeoVolume.h"
55 #include "TVirtualGeoPainter.h"
56 #include "TGeoTube.h"
57 #include "TGeoCone.h"
58 #include "TVirtualPad.h"
59 #include "TBuffer3D.h"
60 #include "TBuffer3DTypes.h"
61 #include "TMath.h"
62 
63 ClassImp(TGeoPcon);
64 
65 ////////////////////////////////////////////////////////////////////////////////
66 /// dummy ctor
67 
68 TGeoPcon::TGeoPcon()
69  :TGeoBBox(),
70  fNz(0),
71  fPhi1(0.),
72  fDphi(0.),
73  fRmin(NULL),
74  fRmax(NULL),
75  fZ(NULL),
76  fFullPhi(kFALSE),
77  fC1(0.),
78  fS1(0.),
79  fC2(0.),
80  fS2(0.),
81  fCm(0.),
82  fSm(0.),
83  fCdphi(0.)
84 {
85  SetShapeBit(TGeoShape::kGeoPcon);
86 }
87 
88 ////////////////////////////////////////////////////////////////////////////////
89 /// Default constructor
90 
91 TGeoPcon::TGeoPcon(Double_t phi, Double_t dphi, Int_t nz)
92  :TGeoBBox(0, 0, 0),
93  fNz(nz),
94  fPhi1(phi),
95  fDphi(dphi),
96  fRmin(NULL),
97  fRmax(NULL),
98  fZ(NULL),
99  fFullPhi(kFALSE),
100  fC1(0.),
101  fS1(0.),
102  fC2(0.),
103  fS2(0.),
104  fCm(0.),
105  fSm(0.),
106  fCdphi(0.)
107 {
108  SetShapeBit(TGeoShape::kGeoPcon);
109  while (fPhi1<0) fPhi1+=360.;
110  fRmin = new Double_t [nz];
111  fRmax = new Double_t [nz];
112  fZ = new Double_t [nz];
113  memset(fRmin, 0, nz*sizeof(Double_t));
114  memset(fRmax, 0, nz*sizeof(Double_t));
115  memset(fZ, 0, nz*sizeof(Double_t));
116  if (TGeoShape::IsSameWithinTolerance(fDphi,360)) fFullPhi = kTRUE;
117  Double_t phi1 = fPhi1;
118  Double_t phi2 = phi1+fDphi;
119  Double_t phim = 0.5*(phi1+phi2);
120  fC1 = TMath::Cos(phi1*TMath::DegToRad());
121  fS1 = TMath::Sin(phi1*TMath::DegToRad());
122  fC2 = TMath::Cos(phi2*TMath::DegToRad());
123  fS2 = TMath::Sin(phi2*TMath::DegToRad());
124  fCm = TMath::Cos(phim*TMath::DegToRad());
125  fSm = TMath::Sin(phim*TMath::DegToRad());
126  fCdphi = TMath::Cos(0.5*fDphi*TMath::DegToRad());
127 }
128 
129 ////////////////////////////////////////////////////////////////////////////////
130 /// Default constructor
131 
132 TGeoPcon::TGeoPcon(const char *name, Double_t phi, Double_t dphi, Int_t nz)
133  :TGeoBBox(name, 0, 0, 0),
134  fNz(nz),
135  fPhi1(phi),
136  fDphi(dphi),
137  fRmin(NULL),
138  fRmax(NULL),
139  fZ(NULL),
140  fFullPhi(kFALSE),
141  fC1(0.),
142  fS1(0.),
143  fC2(0.),
144  fS2(0.),
145  fCm(0.),
146  fSm(0.),
147  fCdphi(0.)
148 {
149  SetShapeBit(TGeoShape::kGeoPcon);
150  while (fPhi1<0) fPhi1+=360.;
151  fRmin = new Double_t [nz];
152  fRmax = new Double_t [nz];
153  fZ = new Double_t [nz];
154  memset(fRmin, 0, nz*sizeof(Double_t));
155  memset(fRmax, 0, nz*sizeof(Double_t));
156  memset(fZ, 0, nz*sizeof(Double_t));
157  if (TGeoShape::IsSameWithinTolerance(fDphi,360)) fFullPhi = kTRUE;
158  Double_t phi1 = fPhi1;
159  Double_t phi2 = phi1+fDphi;
160  Double_t phim = 0.5*(phi1+phi2);
161  fC1 = TMath::Cos(phi1*TMath::DegToRad());
162  fS1 = TMath::Sin(phi1*TMath::DegToRad());
163  fC2 = TMath::Cos(phi2*TMath::DegToRad());
164  fS2 = TMath::Sin(phi2*TMath::DegToRad());
165  fCm = TMath::Cos(phim*TMath::DegToRad());
166  fSm = TMath::Sin(phim*TMath::DegToRad());
167  fCdphi = TMath::Cos(0.5*fDphi*TMath::DegToRad());
168 }
169 
170 ////////////////////////////////////////////////////////////////////////////////
171 /// Default constructor in GEANT3 style
172 /// - param[0] = phi1
173 /// - param[1] = dphi
174 /// - param[2] = nz
175 /// - param[3] = z1
176 /// - param[4] = Rmin1
177 /// - param[5] = Rmax1
178 /// ...
179 
180 TGeoPcon::TGeoPcon(Double_t *param)
181  :TGeoBBox(0, 0, 0),
182  fNz(0),
183  fPhi1(0.),
184  fDphi(0.),
185  fRmin(0),
186  fRmax(0),
187  fZ(0),
188  fFullPhi(kFALSE),
189  fC1(0.),
190  fS1(0.),
191  fC2(0.),
192  fS2(0.),
193  fCm(0.),
194  fSm(0.),
195  fCdphi(0.)
196 {
197  SetShapeBit(TGeoShape::kGeoPcon);
198  SetDimensions(param);
199  ComputeBBox();
200 }
201 
202 ////////////////////////////////////////////////////////////////////////////////
203 ///copy constructor
204 
205 TGeoPcon::TGeoPcon(const TGeoPcon& pc) :
206  TGeoBBox(pc),
207  fNz(0),
208  fPhi1(0.),
209  fDphi(0.),
210  fRmin(0),
211  fRmax(0),
212  fZ(0),
213  fFullPhi(kFALSE),
214  fC1(0.),
215  fS1(0.),
216  fC2(0.),
217  fS2(0.),
218  fCm(0.),
219  fSm(0.),
220  fCdphi(0.)
221 {
222 }
223 
224 ////////////////////////////////////////////////////////////////////////////////
225 ///assignment operator
226 
227 TGeoPcon& TGeoPcon::operator=(const TGeoPcon& pc)
228 {
229  if(this!=&pc) {
230  TGeoBBox::operator=(pc);
231  fNz=0;
232  fPhi1=0.;
233  fDphi=0.;
234  fRmin=0;
235  fRmax=0;
236  fZ=0;
237  fFullPhi=kFALSE;
238  fC1=0;
239  fS1=0;
240  fC2=0;
241  fS2=0;
242  fCm=0;
243  fSm=0;
244  fCdphi=0;
245  }
246  return *this;
247 }
248 
249 ////////////////////////////////////////////////////////////////////////////////
250 /// destructor
251 
252 TGeoPcon::~TGeoPcon()
253 {
254  if (fRmin) {delete[] fRmin; fRmin = 0;}
255  if (fRmax) {delete[] fRmax; fRmax = 0;}
256  if (fZ) {delete[] fZ; fZ = 0;}
257 }
258 
259 ////////////////////////////////////////////////////////////////////////////////
260 /// Computes capacity of the shape in [length^3]
261 
262 Double_t TGeoPcon::Capacity() const
263 {
264  Int_t ipl;
265  Double_t rmin1, rmax1, rmin2, rmax2, phi1, phi2, dz;
266  Double_t capacity = 0.;
267  phi1 = fPhi1;
268  phi2 = fPhi1 + fDphi;
269  for (ipl=0; ipl<fNz-1; ipl++) {
270  dz = 0.5*(fZ[ipl+1]-fZ[ipl]);
271  if (dz < TGeoShape::Tolerance()) continue;
272  rmin1 = fRmin[ipl];
273  rmax1 = fRmax[ipl];
274  rmin2 = fRmin[ipl+1];
275  rmax2 = fRmax[ipl+1];
276  capacity += TGeoConeSeg::Capacity(dz,rmin1,rmax1,rmin2,rmax2,phi1,phi2);
277  }
278  return capacity;
279 }
280 
281 ////////////////////////////////////////////////////////////////////////////////
282 /// compute bounding box of the pcon
283 /// Check if the sections are in increasing Z order
284 
285 void TGeoPcon::ComputeBBox()
286 {
287  for (Int_t isec=0; isec<fNz-1; isec++) {
288  if (TMath::Abs(fZ[isec]-fZ[isec+1]) < TGeoShape::Tolerance()) {
289  fZ[isec+1]=fZ[isec];
290  if (IsSameWithinTolerance(fRmin[isec], fRmin[isec+1]) &&
291  IsSameWithinTolerance(fRmax[isec], fRmax[isec+1])) {
292  InspectShape();
293  Error("ComputeBBox", "Duplicated section %d/%d for shape %s", isec, isec+1, GetName());
294  }
295  }
296  if (fZ[isec]>fZ[isec+1]) {
297  InspectShape();
298  Fatal("ComputeBBox", "Wrong section order");
299  }
300  }
301  // Check if the last sections are valid
302  if (TMath::Abs(fZ[1]-fZ[0]) < TGeoShape::Tolerance() ||
303  TMath::Abs(fZ[fNz-1]-fZ[fNz-2]) < TGeoShape::Tolerance()) {
304  InspectShape();
305  Fatal("ComputeBBox","Shape %s at index %d: Not allowed first two or last two sections at same Z",
306  GetName(), gGeoManager->GetListOfShapes()->IndexOf(this));
307  }
308  Double_t zmin = TMath::Min(fZ[0], fZ[fNz-1]);
309  Double_t zmax = TMath::Max(fZ[0], fZ[fNz-1]);
310  // find largest rmax an smallest rmin
311  Double_t rmin, rmax;
312  rmin = fRmin[TMath::LocMin(fNz, fRmin)];
313  rmax = fRmax[TMath::LocMax(fNz, fRmax)];
314 
315  Double_t xc[4];
316  Double_t yc[4];
317  xc[0] = rmax*fC1;
318  yc[0] = rmax*fS1;
319  xc[1] = rmax*fC2;
320  yc[1] = rmax*fS2;
321  xc[2] = rmin*fC1;
322  yc[2] = rmin*fS1;
323  xc[3] = rmin*fC2;
324  yc[3] = rmin*fS2;
325 
326  Double_t xmin = xc[TMath::LocMin(4, &xc[0])];
327  Double_t xmax = xc[TMath::LocMax(4, &xc[0])];
328  Double_t ymin = yc[TMath::LocMin(4, &yc[0])];
329  Double_t ymax = yc[TMath::LocMax(4, &yc[0])];
330 
331  Double_t ddp = -fPhi1;
332  if (ddp<0) ddp+= 360;
333  if (ddp<=fDphi) xmax = rmax;
334  ddp = 90-fPhi1;
335  if (ddp<0) ddp+= 360;
336  if (ddp<=fDphi) ymax = rmax;
337  ddp = 180-fPhi1;
338  if (ddp<0) ddp+= 360;
339  if (ddp<=fDphi) xmin = -rmax;
340  ddp = 270-fPhi1;
341  if (ddp<0) ddp+= 360;
342  if (ddp<=fDphi) ymin = -rmax;
343  fOrigin[0] = (xmax+xmin)/2;
344  fOrigin[1] = (ymax+ymin)/2;
345  fOrigin[2] = (zmax+zmin)/2;
346  fDX = (xmax-xmin)/2;
347  fDY = (ymax-ymin)/2;
348  fDZ = (zmax-zmin)/2;
349  SetShapeBit(kGeoClosedShape);
350 }
351 
352 ////////////////////////////////////////////////////////////////////////////////
353 /// Compute normal to closest surface from POINT.
354 
355 void TGeoPcon::ComputeNormal(const Double_t *point, const Double_t *dir, Double_t *norm)
356 {
357  memset(norm,0,3*sizeof(Double_t));
358  Double_t r;
359  Double_t ptnew[3];
360  Double_t dz, rmin1, rmax1, rmin2, rmax2;
361  Bool_t is_tube;
362  Int_t ipl = TMath::BinarySearch(fNz, fZ, point[2]);
363  if (ipl==(fNz-1) || ipl<0) {
364  // point outside Z range
365  norm[2] = TMath::Sign(1., dir[2]);
366  return;
367  }
368  Int_t iplclose = ipl;
369  if ((fZ[ipl+1]-point[2])<(point[2]-fZ[ipl])) iplclose++;
370  dz = TMath::Abs(fZ[iplclose]-point[2]);
371  if (dz<1E-5) {
372  if (iplclose==0 || iplclose==(fNz-1)) {
373  norm[2] = TMath::Sign(1., dir[2]);
374  return;
375  }
376  if (iplclose==ipl && TGeoShape::IsSameWithinTolerance(fZ[ipl],fZ[ipl-1])) {
377  r = TMath::Sqrt(point[0]*point[0]+point[1]*point[1]);
378  if (r<TMath::Max(fRmin[ipl],fRmin[ipl-1]) || r>TMath::Min(fRmax[ipl],fRmax[ipl-1])) {
379  norm[2] = TMath::Sign(1., dir[2]);
380  return;
381  }
382  } else {
383  if (TGeoShape::IsSameWithinTolerance(fZ[iplclose],fZ[iplclose+1])) {
384  r = TMath::Sqrt(point[0]*point[0]+point[1]*point[1]);
385  if (r<TMath::Max(fRmin[iplclose],fRmin[iplclose+1]) || r>TMath::Min(fRmax[iplclose],fRmax[iplclose+1])) {
386  norm[2] = TMath::Sign(1., dir[2]);
387  return;
388  }
389  }
390  }
391  } //-> Z done
392  memcpy(ptnew, point, 3*sizeof(Double_t));
393  dz = 0.5*(fZ[ipl+1]-fZ[ipl]);
394  if (TGeoShape::IsSameWithinTolerance(dz,0.)) {
395  norm[2] = TMath::Sign(1., dir[2]);
396  return;
397  }
398  ptnew[2] -= 0.5*(fZ[ipl]+fZ[ipl+1]);
399  rmin1 = fRmin[ipl];
400  rmax1 = fRmax[ipl];
401  rmin2 = fRmin[ipl+1];
402  rmax2 = fRmax[ipl+1];
403  is_tube = (TGeoShape::IsSameWithinTolerance(rmin1,rmin2) && TGeoShape::IsSameWithinTolerance(rmax1,rmax2))?kTRUE:kFALSE;
404  if (!fFullPhi) {
405  if (is_tube) TGeoTubeSeg::ComputeNormalS(ptnew,dir,norm,rmin1,rmax1,dz,fC1,fS1,fC2,fS2);
406  else TGeoConeSeg::ComputeNormalS(ptnew,dir,norm,dz,rmin1,rmax1,rmin2,rmax2,fC1,fS1,fC2,fS2);
407  } else {
408  if (is_tube) TGeoTube::ComputeNormalS(ptnew,dir,norm,rmin1,rmax1,dz);
409  else TGeoCone::ComputeNormalS(ptnew,dir,norm,dz,rmin1,rmax1,rmin2,rmax2);
410  }
411 }
412 
413 ////////////////////////////////////////////////////////////////////////////////
414 /// test if point is inside this shape
415 /// check total z range
416 
417 Bool_t TGeoPcon::Contains(const Double_t *point) const
418 {
419  if ((point[2]<fZ[0]) || (point[2]>fZ[fNz-1])) return kFALSE;
420  // check R squared
421  Double_t r2 = point[0]*point[0]+point[1]*point[1];
422 
423  Int_t izl = 0;
424  Int_t izh = fNz-1;
425  Int_t izt = (fNz-1)/2;
426  while ((izh-izl)>1) {
427  if (point[2] > fZ[izt]) izl = izt;
428  else izh = izt;
429  izt = (izl+izh)>>1;
430  }
431  // the point is in the section bounded by izl and izh Z planes
432 
433  // compute Rmin and Rmax and test the value of R squared
434  Double_t rmin, rmax;
435  if (TGeoShape::IsSameWithinTolerance(fZ[izl],fZ[izh]) && TGeoShape::IsSameWithinTolerance(point[2],fZ[izl])) {
436  rmin = TMath::Min(fRmin[izl], fRmin[izh]);
437  rmax = TMath::Max(fRmax[izl], fRmax[izh]);
438  } else {
439  Double_t dz = fZ[izh] - fZ[izl];
440  Double_t dz1 = point[2] - fZ[izl];
441  rmin = (fRmin[izl]*(dz-dz1)+fRmin[izh]*dz1)/dz;
442  rmax = (fRmax[izl]*(dz-dz1)+fRmax[izh]*dz1)/dz;
443  }
444  if ((r2<rmin*rmin) || (r2>rmax*rmax)) return kFALSE;
445  // now check phi
446  if (TGeoShape::IsSameWithinTolerance(fDphi,360)) return kTRUE;
447  if (r2<1E-10) return kTRUE;
448  Double_t phi = TMath::ATan2(point[1], point[0]) * TMath::RadToDeg();
449  if (phi < 0) phi+=360.0;
450  Double_t ddp = phi-fPhi1;
451  if (ddp<0) ddp+=360.;
452  if (ddp<=fDphi) return kTRUE;
453  return kFALSE;
454 }
455 
456 ////////////////////////////////////////////////////////////////////////////////
457 /// compute closest distance from point px,py to each corner
458 
459 Int_t TGeoPcon::DistancetoPrimitive(Int_t px, Int_t py)
460 {
461  Int_t n = gGeoManager->GetNsegments()+1;
462  const Int_t numPoints = 2*n*fNz;
463  return ShapeDistancetoPrimitive(numPoints, px, py);
464 }
465 
466 ////////////////////////////////////////////////////////////////////////////////
467 /// compute distance from inside point to surface of the polycone
468 
469 Double_t TGeoPcon::DistFromInside(const Double_t *point, const Double_t *dir, Int_t iact, Double_t step, Double_t *safe) const
470 {
471  if (iact<3 && safe) {
472  *safe = Safety(point, kTRUE);
473  if (iact==0) return TGeoShape::Big();
474  if ((iact==1) && (*safe>step)) return TGeoShape::Big();
475  }
476  Double_t snxt = TGeoShape::Big();
477  Double_t sstep = 1E-6;
478  Double_t point_new[3];
479  // determine which z segment contains the point
480  Int_t ipl = TMath::BinarySearch(fNz, fZ, point[2]+TMath::Sign(1.E-10,dir[2]));
481  if (ipl<0) ipl=0;
482  if (ipl==(fNz-1)) ipl--;
483  Double_t dz = 0.5*(fZ[ipl+1]-fZ[ipl]);
484  Bool_t special_case = kFALSE;
485  if (dz<1e-9) {
486  // radius changing segment, make sure track is not in the XY plane
487  if (TGeoShape::IsSameWithinTolerance(dir[2], 0)) {
488  special_case = kTRUE;
489  } else {
490  //check if a close point is still contained
491  point_new[0] = point[0]+sstep*dir[0];
492  point_new[1] = point[1]+sstep*dir[1];
493  point_new[2] = point[2]+sstep*dir[2];
494  if (!Contains(point_new)) return 0.;
495  return (DistFromInside(point_new,dir,iact,step,safe)+sstep);
496  }
497  }
498  // determine if the current segment is a tube or a cone
499  Bool_t intub = kTRUE;
500  if (!TGeoShape::IsSameWithinTolerance(fRmin[ipl],fRmin[ipl+1])) intub=kFALSE;
501  else if (!TGeoShape::IsSameWithinTolerance(fRmax[ipl],fRmax[ipl+1])) intub=kFALSE;
502  // determine phi segmentation
503  memcpy(point_new, point, 2*sizeof(Double_t));
504  // new point in reference system of the current segment
505  point_new[2] = point[2]-0.5*(fZ[ipl]+fZ[ipl+1]);
506 
507  if (special_case) {
508  if (!fFullPhi) snxt = TGeoTubeSeg::DistFromInsideS(point_new, dir,
509  TMath::Min(fRmin[ipl],fRmin[ipl+1]), TMath::Max(fRmax[ipl],fRmax[ipl+1]),
510  dz, fC1,fS1,fC2,fS2,fCm,fSm,fCdphi);
511  else snxt = TGeoTube::DistFromInsideS(point_new, dir,
512  TMath::Min(fRmin[ipl],fRmin[ipl+1]), TMath::Max(fRmax[ipl],fRmax[ipl+1]),dz);
513  return snxt;
514  }
515  if (intub) {
516  if (!fFullPhi) snxt=TGeoTubeSeg::DistFromInsideS(point_new, dir, fRmin[ipl], fRmax[ipl],dz, fC1,fS1,fC2,fS2,fCm,fSm,fCdphi);
517  else snxt=TGeoTube::DistFromInsideS(point_new, dir, fRmin[ipl], fRmax[ipl],dz);
518  } else {
519  if (!fFullPhi) snxt=TGeoConeSeg::DistFromInsideS(point_new,dir,dz,fRmin[ipl],fRmax[ipl],fRmin[ipl+1],fRmax[ipl+1],fC1,fS1,fC2,fS2,fCm,fSm,fCdphi);
520  else snxt=TGeoCone::DistFromInsideS(point_new,dir,dz,fRmin[ipl],fRmax[ipl],fRmin[ipl+1], fRmax[ipl+1]);
521  }
522 
523  for (Int_t i=0; i<3; i++) point_new[i]=point[i]+(snxt+1E-6)*dir[i];
524  if (!Contains(&point_new[0])) return snxt;
525 
526  snxt += DistFromInside(&point_new[0], dir, 3) + 1E-6;
527  return snxt;
528 }
529 
530 ////////////////////////////////////////////////////////////////////////////////
531 /// compute distance to a pcon Z slice. Segment iz must be valid
532 
533 Double_t TGeoPcon::DistToSegZ(const Double_t *point, const Double_t *dir, Int_t &iz) const
534 {
535  Double_t zmin=fZ[iz];
536  Double_t zmax=fZ[iz+1];
537  if (TGeoShape::IsSameWithinTolerance(zmin,zmax)) {
538  if (TGeoShape::IsSameWithinTolerance(dir[2],0)) return TGeoShape::Big();
539  Int_t istep=(dir[2]>0)?1:-1;
540  iz+=istep;
541  if (iz<0 || iz>(fNz-2)) return TGeoShape::Big();
542  return DistToSegZ(point,dir,iz);
543  }
544  Double_t dz=0.5*(zmax-zmin);
545  Double_t local[3];
546  memcpy(&local[0], point, 3*sizeof(Double_t));
547  local[2]=point[2]-0.5*(zmin+zmax);
548  Double_t snxt;
549  Double_t rmin1=fRmin[iz];
550  Double_t rmax1=fRmax[iz];
551  Double_t rmin2=fRmin[iz+1];
552  Double_t rmax2=fRmax[iz+1];
553 
554  if (TGeoShape::IsSameWithinTolerance(rmin1,rmin2) && TGeoShape::IsSameWithinTolerance(rmax1,rmax2)) {
555  if (fFullPhi) snxt=TGeoTube::DistFromOutsideS(local, dir, rmin1, rmax1, dz);
556  else snxt=TGeoTubeSeg::DistFromOutsideS(local,dir,rmin1,rmax1,dz,fC1,fS1,fC2,fS2,fCm,fSm,fCdphi);
557  } else {
558  if (fFullPhi) snxt=TGeoCone::DistFromOutsideS(local,dir,dz,rmin1, rmax1,rmin2,rmax2);
559  else snxt=TGeoConeSeg::DistFromOutsideS(local,dir,dz,rmin1,rmax1,rmin2,rmax2,fC1,fS1,fC2,fS2,fCm,fSm,fCdphi);
560  }
561  if (snxt<1E20) return snxt;
562  // check next segment
563  if (TGeoShape::IsSameWithinTolerance(dir[2],0)) return TGeoShape::Big();
564  Int_t istep=(dir[2]>0)?1:-1;
565  iz+=istep;
566  if (iz<0 || iz>(fNz-2)) return TGeoShape::Big();
567  return DistToSegZ(point,dir,iz);
568 }
569 
570 ////////////////////////////////////////////////////////////////////////////////
571 /// compute distance from outside point to surface of the tube
572 
573 Double_t TGeoPcon::DistFromOutside(const Double_t *point, const Double_t *dir, Int_t iact, Double_t step, Double_t *safe) const
574 {
575  if ((iact<3) && safe) {
576  *safe = Safety(point, kFALSE);
577  if ((iact==1) && (*safe>step)) return TGeoShape::Big();
578  if (iact==0) return TGeoShape::Big();
579  }
580  // check if ray intersect outscribed cylinder
581  if ((point[2]<fZ[0]) && (dir[2]<=0)) return TGeoShape::Big();
582  if ((point[2]>fZ[fNz-1]) && (dir[2]>=0)) return TGeoShape::Big();
583 // Check if the bounding box is crossed within the requested distance
584  Double_t sdist = TGeoBBox::DistFromOutside(point,dir, fDX, fDY, fDZ, fOrigin, step);
585  if (sdist>=step) return TGeoShape::Big();
586 
587  Double_t r2 = point[0]*point[0]+point[1]*point[1];
588  Double_t radmax=0;
589  radmax=fRmax[TMath::LocMax(fNz, fRmax)];
590  if (r2>(radmax*radmax)) {
591  Double_t rpr=-point[0]*dir[0]-point[1]*dir[1];
592  Double_t nxy=dir[0]*dir[0]+dir[1]*dir[1];
593  if (rpr<TMath::Sqrt((r2-radmax*radmax)*nxy)) return TGeoShape::Big();
594  }
595 
596  // find in which Z segment we are
597  Int_t ipl = TMath::BinarySearch(fNz, fZ, point[2]);
598  Int_t ifirst = ipl;
599  if (ifirst<0) {
600  ifirst=0;
601  } else if (ifirst>=(fNz-1)) ifirst=fNz-2;
602  // find if point is in the phi gap
603  Double_t phi=0;
604  if (!fFullPhi) {
605  phi=TMath::ATan2(point[1], point[0]);
606  if (phi<0) phi+=2.*TMath::Pi();
607  }
608 
609  // compute distance to boundary
610  return DistToSegZ(point,dir,ifirst);
611 }
612 
613 ////////////////////////////////////////////////////////////////////////////////
614 /// Defines z position of a section plane, rmin and rmax at this z. Sections
615 /// should be defined in increasing or decreasing Z order and the last section
616 /// HAS to be snum = fNz-1
617 
618 void TGeoPcon::DefineSection(Int_t snum, Double_t z, Double_t rmin, Double_t rmax)
619 {
620  if ((snum<0) || (snum>=fNz)) return;
621  fZ[snum] = z;
622  fRmin[snum] = rmin;
623  fRmax[snum] = rmax;
624  if (rmin>rmax)
625  Warning("DefineSection", "Shape %s: invalid rmin=%g rmax=%g", GetName(), rmin, rmax);
626  if (snum==(fNz-1)) {
627  // Reorder sections in increasing Z order
628  if (fZ[0] > fZ[snum]) {
629  Int_t iz = 0;
630  Int_t izi = fNz-1;
631  Double_t temp;
632  while (iz<izi) {
633  temp = fZ[iz];
634  fZ[iz] = fZ[izi];
635  fZ[izi] = temp;
636  temp = fRmin[iz];
637  fRmin[iz] = fRmin[izi];
638  fRmin[izi] = temp;
639  temp = fRmax[iz];
640  fRmax[iz] = fRmax[izi];
641  fRmax[izi] = temp;
642  iz++;
643  izi--;
644  }
645  }
646  ComputeBBox();
647  }
648 }
649 
650 ////////////////////////////////////////////////////////////////////////////////
651 /// Returns number of segments on each mesh circle segment.
652 
653 Int_t TGeoPcon::GetNsegments() const
654 {
655  return gGeoManager->GetNsegments();
656 }
657 
658 ////////////////////////////////////////////////////////////////////////////////
659 /// Divide this polycone shape belonging to volume "voldiv" into ndiv volumes
660 /// called divname, from start position with the given step. Returns pointer
661 /// to created division cell volume in case of Z divisions. Z divisions can be
662 /// performed if the divided range is in between two consecutive Z planes.
663 /// In case a wrong division axis is supplied, returns pointer to
664 /// volume that was divided.
665 
666 TGeoVolume *TGeoPcon::Divide(TGeoVolume *voldiv, const char *divname, Int_t iaxis, Int_t ndiv,
667  Double_t start, Double_t step)
668 {
669  TGeoShape *shape; //--- shape to be created
670  TGeoVolume *vol; //--- division volume to be created
671  TGeoVolumeMulti *vmulti; //--- generic divided volume
672  TGeoPatternFinder *finder; //--- finder to be attached
673  TString opt = ""; //--- option to be attached
674  Double_t zmin = start;
675  Double_t zmax = start+ndiv*step;
676  Int_t isect = -1;
677  Int_t is, id, ipl;
678  switch (iaxis) {
679  case 1: //--- R division
680  Error("Divide", "Shape %s: cannot divide a pcon on radius", GetName());
681  return 0;
682  case 2: //--- Phi division
683  finder = new TGeoPatternCylPhi(voldiv, ndiv, start, start+ndiv*step);
684  vmulti = gGeoManager->MakeVolumeMulti(divname, voldiv->GetMedium());
685  voldiv->SetFinder(finder);
686  finder->SetDivIndex(voldiv->GetNdaughters());
687  shape = new TGeoPcon(-step/2, step, fNz);
688  for (is=0; is<fNz; is++)
689  ((TGeoPcon*)shape)->DefineSection(is, fZ[is], fRmin[is], fRmax[is]);
690  vol = new TGeoVolume(divname, shape, voldiv->GetMedium());
691  vmulti->AddVolume(vol);
692  opt = "Phi";
693  for (id=0; id<ndiv; id++) {
694  voldiv->AddNodeOffset(vol, id, start+id*step+step/2, opt.Data());
695  ((TGeoNodeOffset*)voldiv->GetNodes()->At(voldiv->GetNdaughters()-1))->SetFinder(finder);
696  }
697  return vmulti;
698  case 3: //--- Z division
699  // find start plane
700  for (ipl=0; ipl<fNz-1; ipl++) {
701  if (start<fZ[ipl]) continue;
702  else {
703  if ((start+ndiv*step)>fZ[ipl+1]) continue;
704  }
705  isect = ipl;
706  zmin = fZ[isect];
707  zmax= fZ[isect+1];
708  break;
709  }
710  if (isect<0) {
711  Error("Divide", "Shape %s: cannot divide pcon on Z if divided region is not between 2 planes", GetName());
712  return 0;
713  }
714  finder = new TGeoPatternZ(voldiv, ndiv, start, start+ndiv*step);
715  vmulti = gGeoManager->MakeVolumeMulti(divname, voldiv->GetMedium());
716  voldiv->SetFinder(finder);
717  finder->SetDivIndex(voldiv->GetNdaughters());
718  opt = "Z";
719  for (id=0; id<ndiv; id++) {
720  Double_t z1 = start+id*step;
721  Double_t z2 = start+(id+1)*step;
722  Double_t rmin1 = (fRmin[isect]*(zmax-z1)-fRmin[isect+1]*(zmin-z1))/(zmax-zmin);
723  Double_t rmax1 = (fRmax[isect]*(zmax-z1)-fRmax[isect+1]*(zmin-z1))/(zmax-zmin);
724  Double_t rmin2 = (fRmin[isect]*(zmax-z2)-fRmin[isect+1]*(zmin-z2))/(zmax-zmin);
725  Double_t rmax2 = (fRmax[isect]*(zmax-z2)-fRmax[isect+1]*(zmin-z2))/(zmax-zmin);
726  Bool_t is_tube = (TGeoShape::IsSameWithinTolerance(fRmin[isect],fRmin[isect+1]) && TGeoShape::IsSameWithinTolerance(fRmax[isect],fRmax[isect+1]))?kTRUE:kFALSE;
727  Bool_t is_seg = (fDphi<360)?kTRUE:kFALSE;
728  if (is_seg) {
729  if (is_tube) shape=new TGeoTubeSeg(fRmin[isect],fRmax[isect],step/2, fPhi1, fPhi1+fDphi);
730  else shape=new TGeoConeSeg(step/2, rmin1, rmax1, rmin2, rmax2, fPhi1, fPhi1+fDphi);
731  } else {
732  if (is_tube) shape=new TGeoTube(fRmin[isect],fRmax[isect],step/2);
733  else shape = new TGeoCone(step/2,rmin1,rmax1,rmin2,rmax2);
734  }
735  vol = new TGeoVolume(divname, shape, voldiv->GetMedium());
736  vmulti->AddVolume(vol);
737  voldiv->AddNodeOffset(vol, id, start+id*step+step/2, opt.Data());
738  ((TGeoNodeOffset*)voldiv->GetNodes()->At(voldiv->GetNdaughters()-1))->SetFinder(finder);
739  }
740  return vmulti;
741  default:
742  Error("Divide", "Shape %s: Wrong axis %d for division",GetName(), iaxis);
743  return 0;
744  }
745 }
746 
747 ////////////////////////////////////////////////////////////////////////////////
748 /// Returns name of axis IAXIS.
749 
750 const char *TGeoPcon::GetAxisName(Int_t iaxis) const
751 {
752  switch (iaxis) {
753  case 1:
754  return "R";
755  case 2:
756  return "PHI";
757  case 3:
758  return "Z";
759  default:
760  return "UNDEFINED";
761  }
762 }
763 
764 ////////////////////////////////////////////////////////////////////////////////
765 /// Get range of shape for a given axis.
766 
767 Double_t TGeoPcon::GetAxisRange(Int_t iaxis, Double_t &xlo, Double_t &xhi) const
768 {
769  xlo = 0;
770  xhi = 0;
771  Double_t dx = 0;
772  switch (iaxis) {
773  case 2:
774  xlo = fPhi1;
775  xhi = fPhi1 + fDphi;
776  dx = fDphi;
777  return dx;
778  case 3:
779  xlo = fZ[0];
780  xhi = fZ[fNz-1];
781  dx = xhi-xlo;
782  return dx;
783  }
784  return dx;
785 }
786 
787 ////////////////////////////////////////////////////////////////////////////////
788 /// Fill vector param[4] with the bounding cylinder parameters. The order
789 /// is the following : Rmin, Rmax, Phi1, Phi2
790 
791 void TGeoPcon::GetBoundingCylinder(Double_t *param) const
792 {
793  param[0] = fRmin[0]; // Rmin
794  param[1] = fRmax[0]; // Rmax
795  for (Int_t i=1; i<fNz; i++) {
796  if (fRmin[i] < param[0]) param[0] = fRmin[i];
797  if (fRmax[i] > param[1]) param[1] = fRmax[i];
798  }
799  param[0] *= param[0];
800  param[1] *= param[1];
801  if (TGeoShape::IsSameWithinTolerance(fDphi,360.)) {
802  param[2] = 0.;
803  param[3] = 360.;
804  return;
805  }
806  param[2] = (fPhi1<0)?(fPhi1+360.):fPhi1; // Phi1
807  param[3] = param[2]+fDphi; // Phi2
808 }
809 
810 ////////////////////////////////////////////////////////////////////////////////
811 /// Returns Rmin for Z segment IPL.
812 
813 Double_t TGeoPcon::GetRmin(Int_t ipl) const
814 {
815  if (ipl<0 || ipl>(fNz-1)) {
816  Error("GetRmin","ipl=%i out of range (0,%i) in shape %s",ipl,fNz-1,GetName());
817  return 0.;
818  }
819  return fRmin[ipl];
820 }
821 
822 ////////////////////////////////////////////////////////////////////////////////
823 /// Returns Rmax for Z segment IPL.
824 
825 Double_t TGeoPcon::GetRmax(Int_t ipl) const
826 {
827  if (ipl<0 || ipl>(fNz-1)) {
828  Error("GetRmax","ipl=%i out of range (0,%i) in shape %s",ipl,fNz-1,GetName());
829  return 0.;
830  }
831  return fRmax[ipl];
832 }
833 
834 ////////////////////////////////////////////////////////////////////////////////
835 /// Returns Z for segment IPL.
836 
837 Double_t TGeoPcon::GetZ(Int_t ipl) const
838 {
839  if (ipl<0 || ipl>(fNz-1)) {
840  Error("GetZ","ipl=%i out of range (0,%i) in shape %s",ipl,fNz-1,GetName());
841  return 0.;
842  }
843  return fZ[ipl];
844 }
845 
846 ////////////////////////////////////////////////////////////////////////////////
847 /// print shape parameters
848 
849 void TGeoPcon::InspectShape() const
850 {
851  printf("*** Shape %s: TGeoPcon ***\n", GetName());
852  printf(" Nz = %i\n", fNz);
853  printf(" phi1 = %11.5f\n", fPhi1);
854  printf(" dphi = %11.5f\n", fDphi);
855  for (Int_t ipl=0; ipl<fNz; ipl++)
856  printf(" plane %i: z=%11.5f Rmin=%11.5f Rmax=%11.5f\n", ipl, fZ[ipl], fRmin[ipl], fRmax[ipl]);
857  printf(" Bounding box:\n");
858  TGeoBBox::InspectShape();
859 }
860 
861 ////////////////////////////////////////////////////////////////////////////////
862 /// Creates a TBuffer3D describing *this* shape.
863 /// Coordinates are in local reference frame.
864 
865 TBuffer3D *TGeoPcon::MakeBuffer3D() const
866 {
867  const Int_t n = gGeoManager->GetNsegments()+1;
868  Int_t nz = GetNz();
869  if (nz < 2) return 0;
870  Int_t nbPnts = nz*2*n;
871  if (nbPnts <= 0) return 0;
872  Double_t dphi = GetDphi();
873 
874  Bool_t specialCase = kFALSE;
875  if (TGeoShape::IsSameWithinTolerance(dphi,360)) specialCase = kTRUE;
876 
877  Int_t nbSegs = 4*(nz*n-1+(specialCase == kTRUE));
878  Int_t nbPols = 2*(nz*n-1+(specialCase == kTRUE));
879  TBuffer3D* buff = new TBuffer3D(TBuffer3DTypes::kGeneric,
880  nbPnts, 3*nbPnts, nbSegs, 3*nbSegs, nbPols, 6*nbPols);
881  if (buff)
882  {
883  SetPoints(buff->fPnts);
884  SetSegsAndPols(*buff);
885  }
886 
887  return buff;
888 }
889 
890 ////////////////////////////////////////////////////////////////////////////////
891 /// Fill TBuffer3D structure for segments and polygons.
892 
893 void TGeoPcon::SetSegsAndPols(TBuffer3D &buff) const
894 {
895  Int_t i, j;
896  const Int_t n = gGeoManager->GetNsegments()+1;
897  Int_t nz = GetNz();
898  if (nz < 2) return;
899  Int_t nbPnts = nz*2*n;
900  if (nbPnts <= 0) return;
901  Double_t dphi = GetDphi();
902 
903  Bool_t specialCase = kFALSE;
904  if (TGeoShape::IsSameWithinTolerance(dphi,360)) specialCase = kTRUE;
905  Int_t c = GetBasicColor();
906 
907  Int_t indx, indx2, k;
908  indx = indx2 = 0;
909 
910  //inside & outside circles, number of segments: 2*nz*(n-1)
911  // special case number of segments: 2*nz*n
912  for (i = 0; i < nz*2; i++) {
913  indx2 = i*n;
914  for (j = 1; j < n; j++) {
915  buff.fSegs[indx++] = c;
916  buff.fSegs[indx++] = indx2+j-1;
917  buff.fSegs[indx++] = indx2+j;
918  }
919  if (specialCase) {
920  buff.fSegs[indx++] = c;
921  buff.fSegs[indx++] = indx2+j-1;
922  buff.fSegs[indx++] = indx2;
923  }
924  }
925 
926  //bottom & top lines, number of segments: 2*n
927  for (i = 0; i < 2; i++) {
928  indx2 = i*(nz-1)*2*n;
929  for (j = 0; j < n; j++) {
930  buff.fSegs[indx++] = c;
931  buff.fSegs[indx++] = indx2+j;
932  buff.fSegs[indx++] = indx2+n+j;
933  }
934  }
935 
936  //inside & outside cilindres, number of segments: 2*(nz-1)*n
937  for (i = 0; i < (nz-1); i++) {
938  //inside cilinder
939  indx2 = i*n*2;
940  for (j = 0; j < n; j++) {
941  buff.fSegs[indx++] = c+2;
942  buff.fSegs[indx++] = indx2+j;
943  buff.fSegs[indx++] = indx2+n*2+j;
944  }
945  //outside cilinder
946  indx2 = i*n*2+n;
947  for (j = 0; j < n; j++) {
948  buff.fSegs[indx++] = c+3;
949  buff.fSegs[indx++] = indx2+j;
950  buff.fSegs[indx++] = indx2+n*2+j;
951  }
952  }
953 
954  //left & right sections, number of segments: 2*(nz-2)
955  // special case number of segments: 0
956  if (!specialCase) {
957  for (i = 1; i < (nz-1); i++) {
958  for (j = 0; j < 2; j++) {
959  buff.fSegs[indx++] = c;
960  buff.fSegs[indx++] = 2*i * n + j*(n-1);
961  buff.fSegs[indx++] = (2*i+1) * n + j*(n-1);
962  }
963  }
964  }
965 
966  Int_t m = n - 1 + (specialCase == kTRUE);
967  indx = 0;
968 
969  //bottom & top, number of polygons: 2*(n-1)
970  // special case number of polygons: 2*n
971  for (j = 0; j < n-1; j++) {
972  buff.fPols[indx++] = c+3;
973  buff.fPols[indx++] = 4;
974  buff.fPols[indx++] = 2*nz*m+j;
975  buff.fPols[indx++] = m+j;
976  buff.fPols[indx++] = 2*nz*m+j+1;
977  buff.fPols[indx++] = j;
978  }
979  for (j = 0; j < n-1; j++) {
980  buff.fPols[indx++] = c+3;
981  buff.fPols[indx++] = 4;
982  buff.fPols[indx++] = 2*nz*m+n+j;
983  buff.fPols[indx++] = (nz*2-2)*m+j;
984  buff.fPols[indx++] = 2*nz*m+n+j+1;
985  buff.fPols[indx++] = (nz*2-2)*m+m+j;
986  }
987  if (specialCase) {
988  buff.fPols[indx++] = c+3;
989  buff.fPols[indx++] = 4;
990  buff.fPols[indx++] = 2*nz*m+j;
991  buff.fPols[indx++] = m+j;
992  buff.fPols[indx++] = 2*nz*m;
993  buff.fPols[indx++] = j;
994 
995  buff.fPols[indx++] = c+3;
996  buff.fPols[indx++] = 4;
997  buff.fPols[indx++] = 2*nz*m+n+j;
998  buff.fPols[indx++] = (nz*2-2)*m+m+j;
999  buff.fPols[indx++] = 2*nz*m+n;
1000  buff.fPols[indx++] = (nz*2-2)*m+j;
1001  }
1002 
1003  //inside & outside, number of polygons: (nz-1)*2*(n-1)
1004  for (k = 0; k < (nz-1); k++) {
1005  for (j = 0; j < n-1; j++) {
1006  buff.fPols[indx++] = c;
1007  buff.fPols[indx++] = 4;
1008  buff.fPols[indx++] = 2*k*m+j;
1009  buff.fPols[indx++] = nz*2*m+(2*k+2)*n+j+1;
1010  buff.fPols[indx++] = (2*k+2)*m+j;
1011  buff.fPols[indx++] = nz*2*m+(2*k+2)*n+j;
1012  }
1013  for (j = 0; j < n-1; j++) {
1014  buff.fPols[indx++] = c+1;
1015  buff.fPols[indx++] = 4;
1016  buff.fPols[indx++] = (2*k+1)*m+j;
1017  buff.fPols[indx++] = nz*2*m+(2*k+3)*n+j;
1018  buff.fPols[indx++] = (2*k+3)*m+j;
1019  buff.fPols[indx++] = nz*2*m+(2*k+3)*n+j+1;
1020  }
1021  if (specialCase) {
1022  buff.fPols[indx++] = c;
1023  buff.fPols[indx++] = 4;
1024  buff.fPols[indx++] = 2*k*m+j;
1025  buff.fPols[indx++] = nz*2*m+(2*k+2)*n;
1026  buff.fPols[indx++] = (2*k+2)*m+j;
1027  buff.fPols[indx++] = nz*2*m+(2*k+2)*n+j;
1028 
1029  buff.fPols[indx++] = c+1;
1030  buff.fPols[indx++] = 4;
1031  buff.fPols[indx++] = (2*k+1)*m+j;
1032  buff.fPols[indx++] = nz*2*m+(2*k+3)*n+j;
1033  buff.fPols[indx++] = (2*k+3)*m+j;
1034  buff.fPols[indx++] = nz*2*m+(2*k+3)*n;
1035  }
1036  }
1037 
1038  //left & right sections, number of polygons: 2*(nz-1)
1039  // special case number of polygons: 0
1040  if (!specialCase) {
1041  indx2 = nz*2*(n-1);
1042  for (k = 0; k < (nz-1); k++) {
1043  buff.fPols[indx++] = c+2;
1044  buff.fPols[indx++] = 4;
1045  buff.fPols[indx++] = k==0 ? indx2 : indx2+2*nz*n+2*(k-1);
1046  buff.fPols[indx++] = indx2+2*(k+1)*n;
1047  buff.fPols[indx++] = indx2+2*nz*n+2*k;
1048  buff.fPols[indx++] = indx2+(2*k+3)*n;
1049 
1050  buff.fPols[indx++] = c+2;
1051  buff.fPols[indx++] = 4;
1052  buff.fPols[indx++] = k==0 ? indx2+n-1 : indx2+2*nz*n+2*(k-1)+1;
1053  buff.fPols[indx++] = indx2+(2*k+3)*n+n-1;
1054  buff.fPols[indx++] = indx2+2*nz*n+2*k+1;
1055  buff.fPols[indx++] = indx2+2*(k+1)*n+n-1;
1056  }
1057  buff.fPols[indx-8] = indx2+n;
1058  buff.fPols[indx-2] = indx2+2*n-1;
1059  }
1060 }
1061 
1062 ////////////////////////////////////////////////////////////////////////////////
1063 /// Compute safety from POINT to segment between planes ipl, ipl+1 within safmin.
1064 
1065 Double_t TGeoPcon::SafetyToSegment(const Double_t *point, Int_t ipl, Bool_t in, Double_t safmin) const
1066 {
1067  Double_t safe = TGeoShape::Big();
1068  if (ipl<0 || ipl>fNz-2) return (safmin+1.); // error in input plane
1069 // Get info about segment.
1070  Double_t dz = 0.5*(fZ[ipl+1]-fZ[ipl]);
1071  if (dz<1E-9) return 1E9; // radius-changing segment
1072  Double_t ptnew[3];
1073  memcpy(ptnew, point, 3*sizeof(Double_t));
1074  ptnew[2] -= 0.5*(fZ[ipl]+fZ[ipl+1]);
1075  safe = TMath::Abs(ptnew[2])-dz;
1076  if (safe>safmin) return TGeoShape::Big(); // means: stop checking further segments
1077  Double_t rmin1 = fRmin[ipl];
1078  Double_t rmax1 = fRmax[ipl];
1079  Double_t rmin2 = fRmin[ipl+1];
1080  Double_t rmax2 = fRmax[ipl+1];
1081  Bool_t is_tube = (TGeoShape::IsSameWithinTolerance(rmin1,rmin2) && TGeoShape::IsSameWithinTolerance(rmax1,rmax2))?kTRUE:kFALSE;
1082  if (!fFullPhi) {
1083  if (is_tube) safe = TGeoTubeSeg::SafetyS(ptnew,in,rmin1,rmax1, dz,fPhi1,fPhi1+fDphi,0);
1084  else safe = TGeoConeSeg::SafetyS(ptnew,in,dz,rmin1,rmax1,rmin2,rmax2,fPhi1,fPhi1+fDphi,0);
1085  } else {
1086  if (is_tube) safe = TGeoTube::SafetyS(ptnew,in,rmin1,rmax1,dz,0);
1087  else safe = TGeoCone::SafetyS(ptnew,in,dz,rmin1,rmax1,rmin2,rmax2,0);
1088  }
1089  if (safe<0) safe=0;
1090  return safe;
1091 }
1092 
1093 ////////////////////////////////////////////////////////////////////////////////
1094 /// computes the closest distance from given point to this shape, according
1095 /// to option. The matching point on the shape is stored in spoint.
1096 /// localize the Z segment
1097 
1098 Double_t TGeoPcon::Safety(const Double_t *point, Bool_t in) const
1099 {
1100  Double_t safmin, saftmp;
1101  Double_t dz;
1102  Int_t ipl, iplane;
1103 
1104  if (in) {
1105  //---> point is inside pcon
1106  ipl = TMath::BinarySearch(fNz, fZ, point[2]);
1107  if (ipl==(fNz-1)) return 0; // point on last Z boundary
1108  if (ipl<0) return 0; // point on first Z boundary
1109  if (ipl>0 && TGeoShape::IsSameWithinTolerance(fZ[ipl-1],fZ[ipl]) && TGeoShape::IsSameWithinTolerance(point[2],fZ[ipl-1])) ipl--;
1110  dz = 0.5*(fZ[ipl+1]-fZ[ipl]);
1111  if (dz<1E-8) {
1112  // Point on a segment-changing plane
1113  safmin = TMath::Min(point[2]-fZ[ipl-1],fZ[ipl+2]-point[2]);
1114  saftmp = TGeoShape::Big();
1115  if (fDphi<360) saftmp = TGeoShape::SafetyPhi(point,in,fPhi1,fPhi1+fDphi);
1116  if (saftmp<safmin) safmin = saftmp;
1117  Double_t radius = TMath::Sqrt(point[0]*point[0]+point[1]*point[1]);
1118  if (fRmin[ipl]>0) safmin = TMath::Min(safmin, radius-fRmin[ipl]);
1119  if (fRmin[ipl+1]>0) safmin = TMath::Min(safmin, radius-fRmin[ipl+1]);
1120  safmin = TMath::Min(safmin, fRmax[ipl]-radius);
1121  safmin = TMath::Min(safmin, fRmax[ipl+1]-radius);
1122  if (safmin<0) safmin = 0;
1123  return safmin;
1124  }
1125  // Check safety for current segment
1126  safmin = SafetyToSegment(point, ipl);
1127  if (safmin>1E10) {
1128  // something went wrong - point is not inside current segment
1129  return 0.;
1130  }
1131  if (safmin<1E-6) return TMath::Abs(safmin); // point on radius-changing plane
1132  // check increasing iplanes
1133  iplane = ipl+1;
1134  saftmp = 0.;
1135 /*
1136  while ((iplane<fNz-1) && saftmp<1E10) {
1137  saftmp = TMath::Abs(SafetyToSegment(point,iplane,kFALSE,safmin));
1138  if (saftmp<safmin) safmin=saftmp;
1139  iplane++;
1140  }
1141  // now decreasing nplanes
1142  iplane = ipl-1;
1143  saftmp = 0.;
1144  while ((iplane>=0) && saftmp<1E10) {
1145  saftmp = TMath::Abs(SafetyToSegment(point,iplane,kFALSE,safmin));
1146  if (saftmp<safmin) safmin=saftmp;
1147  iplane--;
1148  }
1149 */
1150  return safmin;
1151  }
1152  //---> point is outside pcon
1153  ipl = TMath::BinarySearch(fNz, fZ, point[2]);
1154  if (ipl<0) ipl=0;
1155  else if (ipl==fNz-1) ipl=fNz-2;
1156  dz = 0.5*(fZ[ipl+1]-fZ[ipl]);
1157  if (dz<1E-8 && (ipl+2<fNz)) {
1158  ipl++;
1159  dz = 0.5*(fZ[ipl+1]-fZ[ipl]);
1160  }
1161  // Check safety for current segment
1162  safmin = SafetyToSegment(point, ipl, kFALSE);
1163  if (safmin<1E-6) return TMath::Abs(safmin); // point on radius-changing plane
1164  saftmp = 0.;
1165  // check increasing iplanes
1166  iplane = ipl+1;
1167  saftmp = 0.;
1168  while ((iplane<fNz-1) && saftmp<1E10) {
1169  saftmp = TMath::Abs(SafetyToSegment(point,iplane,kFALSE,safmin));
1170  if (saftmp<safmin) safmin=saftmp;
1171  iplane++;
1172  }
1173  // now decreasing nplanes
1174  iplane = ipl-1;
1175  saftmp = 0.;
1176  while ((iplane>=0) && saftmp<1E10) {
1177  saftmp = TMath::Abs(SafetyToSegment(point,iplane,kFALSE,safmin));
1178  if (saftmp<safmin) safmin=saftmp;
1179  iplane--;
1180  }
1181  return safmin;
1182 }
1183 
1184 ////////////////////////////////////////////////////////////////////////////////
1185 /// Save a primitive as a C++ statement(s) on output stream "out".
1186 
1187 void TGeoPcon::SavePrimitive(std::ostream &out, Option_t * /*option*/ /*= ""*/)
1188 {
1189  if (TObject::TestBit(kGeoSavePrimitive)) return;
1190  out << " // Shape: " << GetName() << " type: " << ClassName() << std::endl;
1191  out << " phi1 = " << fPhi1 << ";" << std::endl;
1192  out << " dphi = " << fDphi << ";" << std::endl;
1193  out << " nz = " << fNz << ";" << std::endl;
1194  out << " TGeoPcon *pcon = new TGeoPcon(\"" << GetName() << "\",phi1,dphi,nz);" << std::endl;
1195  for (Int_t i=0; i<fNz; i++) {
1196  out << " z = " << fZ[i] << ";" << std::endl;
1197  out << " rmin = " << fRmin[i] << ";" << std::endl;
1198  out << " rmax = " << fRmax[i] << ";" << std::endl;
1199  out << " pcon->DefineSection(" << i << ", z,rmin,rmax);" << std::endl;
1200  }
1201  out << " TGeoShape *" << GetPointerName() << " = pcon;" << std::endl;
1202  TObject::SetBit(TGeoShape::kGeoSavePrimitive);
1203 }
1204 
1205 ////////////////////////////////////////////////////////////////////////////////
1206 /// Set polycone dimensions starting from an array.
1207 
1208 void TGeoPcon::SetDimensions(Double_t *param)
1209 {
1210  fPhi1 = param[0];
1211  while (fPhi1<0) fPhi1 += 360.;
1212  fDphi = param[1];
1213  fNz = (Int_t)param[2];
1214  if (fNz<2) {
1215  Error("SetDimensions","Pcon %s: Number of Z sections must be > 2", GetName());
1216  return;
1217  }
1218  if (fRmin) delete [] fRmin;
1219  if (fRmax) delete [] fRmax;
1220  if (fZ) delete [] fZ;
1221  fRmin = new Double_t [fNz];
1222  fRmax = new Double_t [fNz];
1223  fZ = new Double_t [fNz];
1224  memset(fRmin, 0, fNz*sizeof(Double_t));
1225  memset(fRmax, 0, fNz*sizeof(Double_t));
1226  memset(fZ, 0, fNz*sizeof(Double_t));
1227  if (TGeoShape::IsSameWithinTolerance(fDphi,360)) fFullPhi = kTRUE;
1228  Double_t phi1 = fPhi1;
1229  Double_t phi2 = phi1+fDphi;
1230  Double_t phim = 0.5*(phi1+phi2);
1231  fC1 = TMath::Cos(phi1*TMath::DegToRad());
1232  fS1 = TMath::Sin(phi1*TMath::DegToRad());
1233  fC2 = TMath::Cos(phi2*TMath::DegToRad());
1234  fS2 = TMath::Sin(phi2*TMath::DegToRad());
1235  fCm = TMath::Cos(phim*TMath::DegToRad());
1236  fSm = TMath::Sin(phim*TMath::DegToRad());
1237  fCdphi = TMath::Cos(0.5*fDphi*TMath::DegToRad());
1238 
1239  for (Int_t i=0; i<fNz; i++)
1240  DefineSection(i, param[3+3*i], param[4+3*i], param[5+3*i]);
1241 }
1242 
1243 ////////////////////////////////////////////////////////////////////////////////
1244 /// create polycone mesh points
1245 
1246 void TGeoPcon::SetPoints(Double_t *points) const
1247 {
1248  Double_t phi, dphi;
1249  Int_t n = gGeoManager->GetNsegments() + 1;
1250  dphi = fDphi/(n-1);
1251  Int_t i, j;
1252  Int_t indx = 0;
1253 
1254  if (points) {
1255  for (i = 0; i < fNz; i++) {
1256  for (j = 0; j < n; j++) {
1257  phi = (fPhi1+j*dphi)*TMath::DegToRad();
1258  points[indx++] = fRmin[i] * TMath::Cos(phi);
1259  points[indx++] = fRmin[i] * TMath::Sin(phi);
1260  points[indx++] = fZ[i];
1261  }
1262  for (j = 0; j < n; j++) {
1263  phi = (fPhi1+j*dphi)*TMath::DegToRad();
1264  points[indx++] = fRmax[i] * TMath::Cos(phi);
1265  points[indx++] = fRmax[i] * TMath::Sin(phi);
1266  points[indx++] = fZ[i];
1267  }
1268  }
1269  }
1270 }
1271 
1272 ////////////////////////////////////////////////////////////////////////////////
1273 /// create polycone mesh points
1274 
1275 void TGeoPcon::SetPoints(Float_t *points) const
1276 {
1277  Double_t phi, dphi;
1278  Int_t n = gGeoManager->GetNsegments() + 1;
1279  dphi = fDphi/(n-1);
1280  Int_t i, j;
1281  Int_t indx = 0;
1282 
1283  if (points) {
1284  for (i = 0; i < fNz; i++) {
1285  for (j = 0; j < n; j++) {
1286  phi = (fPhi1+j*dphi)*TMath::DegToRad();
1287  points[indx++] = fRmin[i] * TMath::Cos(phi);
1288  points[indx++] = fRmin[i] * TMath::Sin(phi);
1289  points[indx++] = fZ[i];
1290  }
1291  for (j = 0; j < n; j++) {
1292  phi = (fPhi1+j*dphi)*TMath::DegToRad();
1293  points[indx++] = fRmax[i] * TMath::Cos(phi);
1294  points[indx++] = fRmax[i] * TMath::Sin(phi);
1295  points[indx++] = fZ[i];
1296  }
1297  }
1298  }
1299 }
1300 ////////////////////////////////////////////////////////////////////////////////
1301 /// Return number of vertices of the mesh representation
1302 
1303 Int_t TGeoPcon::GetNmeshVertices() const
1304 {
1305  Int_t n = gGeoManager->GetNsegments()+1;
1306  Int_t numPoints = fNz*2*n;
1307  return numPoints;
1308 }
1309 
1310 ////////////////////////////////////////////////////////////////////////////////
1311 /// fill size of this 3-D object
1312 
1313 void TGeoPcon::Sizeof3D() const
1314 {
1315 }
1316 
1317 ////////////////////////////////////////////////////////////////////////////////
1318 /// Returns numbers of vertices, segments and polygons composing the shape mesh.
1319 
1320 void TGeoPcon::GetMeshNumbers(Int_t &nvert, Int_t &nsegs, Int_t &npols) const
1321 {
1322  Int_t n = gGeoManager->GetNsegments()+1;
1323  Int_t nz = GetNz();
1324  nvert = nz*2*n;
1325  Bool_t specialCase = TGeoShape::IsSameWithinTolerance(GetDphi(),360);
1326  nsegs = 4*(nz*n-1+(specialCase == kTRUE));
1327  npols = 2*(nz*n-1+(specialCase == kTRUE));
1328 }
1329 
1330 ////////////////////////////////////////////////////////////////////////////////
1331 /// Fills a static 3D buffer and returns a reference.
1332 
1333 const TBuffer3D & TGeoPcon::GetBuffer3D(Int_t reqSections, Bool_t localFrame) const
1334 {
1335  static TBuffer3D buffer(TBuffer3DTypes::kGeneric);
1336 
1337  TGeoBBox::FillBuffer3D(buffer, reqSections, localFrame);
1338 
1339  if (reqSections & TBuffer3D::kRawSizes) {
1340  const Int_t n = gGeoManager->GetNsegments()+1;
1341  Int_t nz = GetNz();
1342  Int_t nbPnts = nz*2*n;
1343  if (nz >= 2 && nbPnts > 0) {
1344  Bool_t specialCase = TGeoShape::IsSameWithinTolerance(GetDphi(),360);
1345  Int_t nbSegs = 4*(nz*n-1+(specialCase == kTRUE));
1346  Int_t nbPols = 2*(nz*n-1+(specialCase == kTRUE));
1347  if (buffer.SetRawSizes(nbPnts, 3*nbPnts, nbSegs, 3*nbSegs, nbPols, 6*nbPols)) {
1348  buffer.SetSectionsValid(TBuffer3D::kRawSizes);
1349  }
1350  }
1351  }
1352  // TODO: Push down to TGeoShape?? Would have to do raw sizes set first..
1353  // can rest of TGeoShape be deferred until after this?
1354  if ((reqSections & TBuffer3D::kRaw) && buffer.SectionsValid(TBuffer3D::kRawSizes)) {
1355  SetPoints(buffer.fPnts);
1356  if (!buffer.fLocalFrame) {
1357  TransformPoints(buffer.fPnts, buffer.NbPnts());
1358  }
1359 
1360  SetSegsAndPols(buffer);
1361  buffer.SetSectionsValid(TBuffer3D::kRaw);
1362  }
1363 
1364  return buffer;
1365 }
1366 
1367 ////////////////////////////////////////////////////////////////////////////////
1368 /// Stream an object of class TGeoPcon.
1369 
1370 void TGeoPcon::Streamer(TBuffer &R__b)
1371 {
1372  if (R__b.IsReading()) {
1373  R__b.ReadClassBuffer(TGeoPcon::Class(),this);
1374  if (TGeoShape::IsSameWithinTolerance(fDphi,360)) fFullPhi = kTRUE;
1375  Double_t phi1 = fPhi1;
1376  Double_t phi2 = phi1+fDphi;
1377  Double_t phim = 0.5*(phi1+phi2);
1378  fC1 = TMath::Cos(phi1*TMath::DegToRad());
1379  fS1 = TMath::Sin(phi1*TMath::DegToRad());
1380  fC2 = TMath::Cos(phi2*TMath::DegToRad());
1381  fS2 = TMath::Sin(phi2*TMath::DegToRad());
1382  fCm = TMath::Cos(phim*TMath::DegToRad());
1383  fSm = TMath::Sin(phim*TMath::DegToRad());
1384  fCdphi = TMath::Cos(0.5*fDphi*TMath::DegToRad());
1385  } else {
1386  R__b.WriteClassBuffer(TGeoPcon::Class(),this);
1387  }
1388 }
1389 
1390 ////////////////////////////////////////////////////////////////////////////////
1391 /// Check the inside status for each of the points in the array.
1392 /// Input: Array of point coordinates + vector size
1393 /// Output: Array of Booleans for the inside of each point
1394 
1395 void TGeoPcon::Contains_v(const Double_t *points, Bool_t *inside, Int_t vecsize) const
1396 {
1397  for (Int_t i=0; i<vecsize; i++) inside[i] = Contains(&points[3*i]);
1398 }
1399 
1400 ////////////////////////////////////////////////////////////////////////////////
1401 /// Compute the normal for an array o points so that norm.dot.dir is positive
1402 /// Input: Arrays of point coordinates and directions + vector size
1403 /// Output: Array of normal directions
1404 
1405 void TGeoPcon::ComputeNormal_v(const Double_t *points, const Double_t *dirs, Double_t *norms, Int_t vecsize)
1406 {
1407  for (Int_t i=0; i<vecsize; i++) ComputeNormal(&points[3*i], &dirs[3*i], &norms[3*i]);
1408 }
1409 
1410 ////////////////////////////////////////////////////////////////////////////////
1411 /// Compute distance from array of input points having directions specified by dirs. Store output in dists
1412 
1413 void TGeoPcon::DistFromInside_v(const Double_t *points, const Double_t *dirs, Double_t *dists, Int_t vecsize, Double_t* step) const
1414 {
1415  for (Int_t i=0; i<vecsize; i++) dists[i] = DistFromInside(&points[3*i], &dirs[3*i], 3, step[i]);
1416 }
1417 
1418 ////////////////////////////////////////////////////////////////////////////////
1419 /// Compute distance from array of input points having directions specified by dirs. Store output in dists
1420 
1421 void TGeoPcon::DistFromOutside_v(const Double_t *points, const Double_t *dirs, Double_t *dists, Int_t vecsize, Double_t* step) const
1422 {
1423  for (Int_t i=0; i<vecsize; i++) dists[i] = DistFromOutside(&points[3*i], &dirs[3*i], 3, step[i]);
1424 }
1425 
1426 ////////////////////////////////////////////////////////////////////////////////
1427 /// Compute safe distance from each of the points in the input array.
1428 /// Input: Array of point coordinates, array of statuses for these points, size of the arrays
1429 /// Output: Safety values
1430 
1431 void TGeoPcon::Safety_v(const Double_t *points, const Bool_t *inside, Double_t *safe, Int_t vecsize) const
1432 {
1433  for (Int_t i=0; i<vecsize; i++) safe[i] = Safety(&points[3*i], inside[i]);
1434 }