Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
TGeoPgon.cxx
Go to the documentation of this file.
1 // @(#)root/geom:$Id$
2 // Author: Andrei Gheata 31/01/02
3 // TGeoPgon::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 TGeoPgon
14 \ingroup Geometry_classes
15 
16 A polygone. It has at least 10 parameters :
17  - the lower phi limit;
18  - the range in phi;
19  - the number of equal edges on each z plane;
20  - the number of z planes (at least two) where the inner/outer
21  radii are changing;
22  - z coordinate, inner and outer radius of the inscribed cercle
23  (distance from center to edges) for each z plane
24 
25 Begin_Macro(source)
26 {
27  TCanvas *c = new TCanvas("c", "c",0,0,600,600);
28  new TGeoManager("pgon", "poza11");
29  TGeoMaterial *mat = new TGeoMaterial("Al", 26.98,13,2.7);
30  TGeoMedium *med = new TGeoMedium("MED",1,mat);
31  TGeoVolume *top = gGeoManager->MakeBox("TOP",med,150,150,100);
32  gGeoManager->SetTopVolume(top);
33  TGeoVolume *vol = gGeoManager->MakePgon("PGON",med, -45.0,270.0,4,4);
34  TGeoPgon *pgon = (TGeoPgon*)(vol->GetShape());
35  pgon->DefineSection(0,-70,45,50);
36  pgon->DefineSection(1,0,35,40);
37  pgon->DefineSection(2,0,30,35);
38  pgon->DefineSection(3,70,90,100);
39  vol->SetLineWidth(2);
40  top->AddNode(vol,1);
41  gGeoManager->CloseGeometry();
42  gGeoManager->SetNsegments(80);
43  top->Draw();
44  TView *view = gPad->GetView();
45  view->ShowAxis();
46 }
47 End_Macro
48 */
49 
50 #include "Riostream.h"
51 
52 #include "TGeoPgon.h"
53 
54 #include "TGeoManager.h"
55 #include "TGeoVolume.h"
56 #include "TVirtualGeoPainter.h"
57 #include "TGeoTube.h"
58 #include "TVirtualPad.h"
59 #include "TBuffer3D.h"
60 #include "TBuffer3DTypes.h"
61 #include "TMath.h"
62 
63 ClassImp(TGeoPgon);
64 
65  ////////////////////////////////////////////////////////////////////////////////
66  /// Constructor.
67 
68  TGeoPgon::ThreadData_t::ThreadData_t()
69  : fIntBuffer(0), fDblBuffer(0)
70 {
71 }
72 
73 ////////////////////////////////////////////////////////////////////////////////
74 /// Destructor.
75 
76 TGeoPgon::ThreadData_t::~ThreadData_t()
77 {
78  delete[] fIntBuffer;
79  delete[] fDblBuffer;
80 }
81 
82 ////////////////////////////////////////////////////////////////////////////////
83 
84 TGeoPgon::ThreadData_t &TGeoPgon::GetThreadData() const
85 {
86  Int_t tid = TGeoManager::ThreadId();
87  return *fThreadData[tid];
88 }
89 
90 ////////////////////////////////////////////////////////////////////////////////
91 
92 void TGeoPgon::ClearThreadData() const
93 {
94  std::lock_guard<std::mutex> guard(fMutex);
95  std::vector<ThreadData_t *>::iterator i = fThreadData.begin();
96  while (i != fThreadData.end()) {
97  delete *i;
98  ++i;
99  }
100  fThreadData.clear();
101  fThreadSize = 0;
102 }
103 
104 ////////////////////////////////////////////////////////////////////////////////
105 /// Create thread data for n threads max.
106 
107 void TGeoPgon::CreateThreadData(Int_t nthreads)
108 {
109  if (fThreadSize) ClearThreadData();
110  std::lock_guard<std::mutex> guard(fMutex);
111  fThreadData.resize(nthreads);
112  fThreadSize = nthreads;
113  for (Int_t tid = 0; tid < nthreads; tid++) {
114  if (fThreadData[tid] == 0) {
115  fThreadData[tid] = new ThreadData_t;
116  fThreadData[tid]->fIntBuffer = new Int_t[fNedges + 10];
117  fThreadData[tid]->fDblBuffer = new Double_t[fNedges + 10];
118  }
119  }
120 }
121 
122 ////////////////////////////////////////////////////////////////////////////////
123 /// dummy ctor
124 
125 TGeoPgon::TGeoPgon()
126 {
127  SetShapeBit(TGeoShape::kGeoPgon);
128  fNedges = 0;
129  fThreadSize = 0;
130 }
131 
132 ////////////////////////////////////////////////////////////////////////////////
133 /// Default constructor
134 
135 TGeoPgon::TGeoPgon(Double_t phi, Double_t dphi, Int_t nedges, Int_t nz) : TGeoPcon(phi, dphi, nz)
136 {
137  SetShapeBit(TGeoShape::kGeoPgon);
138  fNedges = nedges;
139  fThreadSize = 0;
140  CreateThreadData(1);
141 }
142 
143 ////////////////////////////////////////////////////////////////////////////////
144 /// Default constructor
145 
146 TGeoPgon::TGeoPgon(const char *name, Double_t phi, Double_t dphi, Int_t nedges, Int_t nz)
147  : TGeoPcon(name, phi, dphi, nz)
148 {
149  SetShapeBit(TGeoShape::kGeoPgon);
150  fNedges = nedges;
151  fThreadSize = 0;
152  CreateThreadData(1);
153 }
154 
155 ////////////////////////////////////////////////////////////////////////////////
156 /// Default constructor in GEANT3 style
157 /// - param[0] = phi1
158 /// - param[1] = dphi
159 /// - param[2] = nedges
160 /// - param[3] = nz
161 /// - param[4] = z1
162 /// - param[5] = Rmin1
163 /// - param[6] = Rmax1
164 /// ...
165 
166 TGeoPgon::TGeoPgon(Double_t *param) : TGeoPcon()
167 {
168  SetShapeBit(TGeoShape::kGeoPgon);
169  SetDimensions(param);
170  ComputeBBox();
171  fThreadSize = 0;
172  CreateThreadData(1);
173 }
174 
175 ////////////////////////////////////////////////////////////////////////////////
176 /// destructor
177 
178 TGeoPgon::~TGeoPgon()
179 {
180  ClearThreadData();
181 }
182 
183 ////////////////////////////////////////////////////////////////////////////////
184 /// Computes capacity of the shape in [length^3]
185 
186 Double_t TGeoPgon::Capacity() const
187 {
188  Int_t ipl;
189  Double_t rmin1, rmax1, rmin2, rmax2, dphi, dz;
190  Double_t capacity = 0.;
191  dphi = fDphi / fNedges; // [deg]
192  Double_t tphi2 = TMath::Tan(0.5 * dphi * TMath::DegToRad());
193  for (ipl = 0; ipl < fNz - 1; ipl++) {
194  dz = fZ[ipl + 1] - fZ[ipl];
195  if (dz < TGeoShape::Tolerance()) continue;
196  rmin1 = fRmin[ipl];
197  rmax1 = fRmax[ipl];
198  rmin2 = fRmin[ipl + 1];
199  rmax2 = fRmax[ipl + 1];
200  capacity += fNedges * (tphi2 / 3.) * dz *
201  (rmax1 * rmax1 + rmax1 * rmax2 + rmax2 * rmax2 - rmin1 * rmin1 - rmin1 * rmin2 - rmin2 * rmin2);
202  }
203  return capacity;
204 }
205 
206 ////////////////////////////////////////////////////////////////////////////////
207 /// compute bounding box for a polygone
208 /// Check if the sections are in increasing Z order
209 
210 void TGeoPgon::ComputeBBox()
211 {
212  for (Int_t isec = 0; isec < fNz - 1; isec++) {
213  if (fZ[isec] > fZ[isec + 1]) {
214  InspectShape();
215  Fatal("ComputeBBox", "Wrong section order");
216  }
217  }
218  // Check if the last sections are valid
219  if (TMath::Abs(fZ[1] - fZ[0]) < TGeoShape::Tolerance() ||
220  TMath::Abs(fZ[fNz - 1] - fZ[fNz - 2]) < TGeoShape::Tolerance()) {
221  InspectShape();
222  Fatal("ComputeBBox", "Shape %s at index %d: Not allowed first two or last two sections at same Z", GetName(),
223  gGeoManager->GetListOfShapes()->IndexOf(this));
224  }
225  Double_t zmin = TMath::Min(fZ[0], fZ[fNz - 1]);
226  Double_t zmax = TMath::Max(fZ[0], fZ[fNz - 1]);
227  // find largest rmax an smallest rmin
228  Double_t rmin, rmax;
229  Double_t divphi = fDphi / fNedges;
230  // find the radius of the outscribed circle
231  rmin = fRmin[TMath::LocMin(fNz, fRmin)];
232  rmax = fRmax[TMath::LocMax(fNz, fRmax)];
233  rmax = rmax / TMath::Cos(0.5 * divphi * TMath::DegToRad());
234  Double_t phi1 = fPhi1;
235  Double_t phi2 = phi1 + fDphi;
236 
237  Double_t xc[4];
238  Double_t yc[4];
239  xc[0] = rmax * TMath::Cos(phi1 * TMath::DegToRad());
240  yc[0] = rmax * TMath::Sin(phi1 * TMath::DegToRad());
241  xc[1] = rmax * TMath::Cos(phi2 * TMath::DegToRad());
242  yc[1] = rmax * TMath::Sin(phi2 * TMath::DegToRad());
243  xc[2] = rmin * TMath::Cos(phi1 * TMath::DegToRad());
244  yc[2] = rmin * TMath::Sin(phi1 * TMath::DegToRad());
245  xc[3] = rmin * TMath::Cos(phi2 * TMath::DegToRad());
246  yc[3] = rmin * TMath::Sin(phi2 * TMath::DegToRad());
247 
248  Double_t xmin = xc[TMath::LocMin(4, &xc[0])];
249  Double_t xmax = xc[TMath::LocMax(4, &xc[0])];
250  Double_t ymin = yc[TMath::LocMin(4, &yc[0])];
251  Double_t ymax = yc[TMath::LocMax(4, &yc[0])];
252 
253  Double_t ddp = -phi1;
254  if (ddp < 0) ddp += 360;
255  if (ddp <= fDphi) xmax = rmax;
256  ddp = 90 - phi1;
257  if (ddp < 0) ddp += 360;
258  if (ddp <= fDphi) ymax = rmax;
259  ddp = 180 - phi1;
260  if (ddp < 0) ddp += 360;
261  if (ddp <= fDphi) xmin = -rmax;
262  ddp = 270 - phi1;
263  if (ddp < 0) ddp += 360;
264  if (ddp <= fDphi) ymin = -rmax;
265  fOrigin[0] = 0.5 * (xmax + xmin);
266  fOrigin[1] = 0.5 * (ymax + ymin);
267  fOrigin[2] = 0.5 * (zmax + zmin);
268  fDX = 0.5 * (xmax - xmin);
269  fDY = 0.5 * (ymax - ymin);
270  fDZ = 0.5 * (zmax - zmin);
271  SetShapeBit(kGeoClosedShape);
272 }
273 
274 ////////////////////////////////////////////////////////////////////////////////
275 /// Compute normal to closest surface from POINT.
276 
277 void TGeoPgon::ComputeNormal(const Double_t *point, const Double_t *dir, Double_t *norm)
278 {
279  memset(norm, 0, 3 * sizeof(Double_t));
280  Double_t phi1 = 0, phi2 = 0, c1 = 0, s1 = 0, c2 = 0, s2 = 0;
281  Double_t dz, rmin1, rmin2;
282  Bool_t is_seg = (fDphi < 360) ? kTRUE : kFALSE;
283  if (is_seg) {
284  phi1 = fPhi1;
285  if (phi1 < 0) phi1 += 360;
286  phi2 = phi1 + fDphi;
287  phi1 *= TMath::DegToRad();
288  phi2 *= TMath::DegToRad();
289  c1 = TMath::Cos(phi1);
290  s1 = TMath::Sin(phi1);
291  c2 = TMath::Cos(phi2);
292  s2 = TMath::Sin(phi2);
293  if (TGeoShape::IsCloseToPhi(1E-5, point, c1, s1, c2, s2)) {
294  TGeoShape::NormalPhi(point, dir, norm, c1, s1, c2, s2);
295  return;
296  }
297  } // Phi done
298 
299  Int_t ipl = TMath::BinarySearch(fNz, fZ, point[2]);
300  if (ipl == (fNz - 1) || ipl < 0) {
301  // point outside Z range
302  norm[2] = TMath::Sign(1., dir[2]);
303  return;
304  }
305  Int_t iplclose = ipl;
306  if ((fZ[ipl + 1] - point[2]) < (point[2] - fZ[ipl])) iplclose++;
307  dz = TMath::Abs(fZ[iplclose] - point[2]);
308 
309  Double_t divphi = fDphi / fNedges;
310  Double_t phi = TMath::ATan2(point[1], point[0]) * TMath::RadToDeg();
311  while (phi < fPhi1) phi += 360.;
312  Double_t ddp = phi - fPhi1;
313  Int_t ipsec = Int_t(ddp / divphi);
314  Double_t ph0 = (fPhi1 + divphi * (ipsec + 0.5)) * TMath::DegToRad();
315  // compute projected distance
316  Double_t r, rsum, rpgon, ta, calf;
317  r = TMath::Abs(point[0] * TMath::Cos(ph0) + point[1] * TMath::Sin(ph0));
318  if (dz < 1E-5) {
319  if (iplclose == 0 || iplclose == (fNz - 1)) {
320  norm[2] = TMath::Sign(1., dir[2]);
321  return;
322  }
323  if (iplclose == ipl && TGeoShape::IsSameWithinTolerance(fZ[ipl], fZ[ipl - 1])) {
324  if (r < TMath::Max(fRmin[ipl], fRmin[ipl - 1]) || r > TMath::Min(fRmax[ipl], fRmax[ipl - 1])) {
325  norm[2] = TMath::Sign(1., dir[2]);
326  return;
327  }
328  } else {
329  if (TGeoShape::IsSameWithinTolerance(fZ[iplclose], fZ[iplclose + 1])) {
330  if (r < TMath::Max(fRmin[iplclose], fRmin[iplclose + 1]) ||
331  r > TMath::Min(fRmax[iplclose], fRmax[iplclose + 1])) {
332  norm[2] = TMath::Sign(1., dir[2]);
333  return;
334  }
335  }
336  }
337  } //-> Z done
338 
339  dz = fZ[ipl + 1] - fZ[ipl];
340  rmin1 = fRmin[ipl];
341  rmin2 = fRmin[ipl + 1];
342  rsum = rmin1 + rmin2;
343  Double_t safe = TGeoShape::Big();
344  if (rsum > 1E-10) {
345  ta = (rmin2 - rmin1) / dz;
346  calf = 1. / TMath::Sqrt(1 + ta * ta);
347  rpgon = rmin1 + (point[2] - fZ[ipl]) * ta;
348  safe = TMath::Abs(r - rpgon);
349  norm[0] = calf * TMath::Cos(ph0);
350  norm[1] = calf * TMath::Sin(ph0);
351  norm[2] = -calf * ta;
352  }
353  ta = (fRmax[ipl + 1] - fRmax[ipl]) / dz;
354  calf = 1. / TMath::Sqrt(1 + ta * ta);
355  rpgon = fRmax[ipl] + (point[2] - fZ[ipl]) * ta;
356  if (safe > TMath::Abs(rpgon - r)) {
357  norm[0] = calf * TMath::Cos(ph0);
358  norm[1] = calf * TMath::Sin(ph0);
359  norm[2] = -calf * ta;
360  }
361  if (norm[0] * dir[0] + norm[1] * dir[1] + norm[2] * dir[2] < 0) {
362  norm[0] = -norm[0];
363  norm[1] = -norm[1];
364  norm[2] = -norm[2];
365  }
366 }
367 
368 ////////////////////////////////////////////////////////////////////////////////
369 /// test if point is inside this shape
370 /// check total z range
371 
372 Bool_t TGeoPgon::Contains(const Double_t *point) const
373 {
374  if (point[2] < fZ[0]) return kFALSE;
375  if (point[2] > fZ[fNz - 1]) return kFALSE;
376  Double_t divphi = fDphi / fNedges;
377  // now check phi
378  Double_t phi = TMath::ATan2(point[1], point[0]) * TMath::RadToDeg();
379  while (phi < fPhi1) phi += 360.0;
380  Double_t ddp = phi - fPhi1;
381  if (ddp > fDphi) return kFALSE;
382  // now find phi division
383  Int_t ipsec = TMath::Min(Int_t(ddp / divphi), fNedges - 1);
384  Double_t ph0 = (fPhi1 + divphi * (ipsec + 0.5)) * TMath::DegToRad();
385  // now check projected distance
386  Double_t r = point[0] * TMath::Cos(ph0) + point[1] * TMath::Sin(ph0);
387  // find in which Z section the point is in
388  Int_t iz = TMath::BinarySearch(fNz, fZ, point[2]);
389  if (iz == fNz - 1) {
390  if (r < fRmin[iz]) return kFALSE;
391  if (r > fRmax[iz]) return kFALSE;
392  return kTRUE;
393  }
394  Double_t dz = fZ[iz + 1] - fZ[iz];
395  Double_t rmin, rmax;
396  if (dz < 1E-8) {
397  // we are at a radius-changing plane
398  rmin = TMath::Min(fRmin[iz], fRmin[iz + 1]);
399  rmax = TMath::Max(fRmax[iz], fRmax[iz + 1]);
400  if (r < rmin) return kFALSE;
401  if (r > rmax) return kFALSE;
402  return kTRUE;
403  }
404  // now compute rmin and rmax and test the value of r
405  Double_t dzrat = (point[2] - fZ[iz]) / dz;
406  rmin = fRmin[iz] + dzrat * (fRmin[iz + 1] - fRmin[iz]);
407  // is the point inside the 'hole' at the center of the volume ?
408  if (r < rmin) return kFALSE;
409  rmax = fRmax[iz] + dzrat * (fRmax[iz + 1] - fRmax[iz]);
410  if (r > rmax) return kFALSE;
411 
412  return kTRUE;
413 }
414 
415 ////////////////////////////////////////////////////////////////////////////////
416 /// compute distance from inside point to surface of the polygone
417 /// first find out in which Z section the point is in
418 
419 Double_t TGeoPgon::DistFromInside(const Double_t *point, const Double_t *dir, Int_t iact, Double_t step,
420  Double_t *safe) const
421 {
422  if (iact < 3 && safe) {
423  *safe = Safety(point, kTRUE);
424  if (iact == 0) return TGeoShape::Big();
425  if (iact == 1 && step < *safe) return TGeoShape::Big();
426  }
427  // find current Z section
428  Int_t ipl, ipsec;
429  ipl = TMath::BinarySearch(fNz, fZ, point[2]);
430  if (ipl == fNz - 1) {
431  if (dir[2] >= 0) return 0.;
432  ipl--;
433  }
434  if (ipl < 0) {
435  // point out
436  if (dir[2] <= 0) return 0.;
437  ipl++;
438  }
439  Double_t stepmax = step;
440  if (!fThreadSize) ((TGeoPgon *)this)->CreateThreadData(1);
441  ThreadData_t &td = GetThreadData();
442  Double_t *sph = td.fDblBuffer;
443  Int_t *iph = td.fIntBuffer;
444  // locate current phi sector [0,fNedges-1]; -1 for dead region
445  LocatePhi(point, ipsec);
446  if (ipsec < 0) {
447  // Point on a phi boundary - entering or exiting ?
448  Double_t phi1 = fPhi1 * TMath::DegToRad();
449  Double_t phi2 = (fPhi1 + fDphi) * TMath::DegToRad();
450  if ((point[0] * dir[1] - point[1] * dir[0]) > 0) {
451  // phi1 next crossing
452  if ((point[0] * TMath::Cos(phi1) + point[1] * TMath::Sin(phi1)) <
453  (point[0] * TMath::Cos(phi2) + point[1] * TMath::Sin(phi2))) {
454  // close to phimax
455  return 0.0;
456  } else {
457  // close to phi1 - ignore it
458  ipsec = 0;
459  }
460  } else {
461  // phimax next crossing
462  if ((point[0] * TMath::Cos(phi1) + point[1] * TMath::Sin(phi1)) >
463  (point[0] * TMath::Cos(phi2) + point[1] * TMath::Sin(phi2))) {
464  // close to phi1
465  return 0.0;
466  } else {
467  // close to phimax - ignore it
468  ipsec = fNedges - 1;
469  }
470  }
471  }
472  Int_t ipln = -1;
473  if (TGeoShape::IsSameWithinTolerance(fZ[ipl], fZ[ipl + 1])) {
474  ipln = ipl;
475  } else {
476  if (fNz > 3 && ipl >= 0 && ipl < fNz - 3 && TGeoShape::IsSameWithinTolerance(fZ[ipl + 1], fZ[ipl + 2]) &&
477  TMath::Abs(point[2] - fZ[ipl + 1]) < 1.E-8) {
478  ipln = ipl + 1;
479  } else {
480  if (ipl > 1 && TGeoShape::IsSameWithinTolerance(fZ[ipl], fZ[ipl - 1]) &&
481  TMath::Abs(point[2] - fZ[ipl]) < 1.E-8)
482  ipln = ipl - 1;
483  }
484  }
485  if (ipln > 0) {
486  // point between segments
487  Double_t divphi = fDphi / fNedges;
488  Double_t phi = (fPhi1 + (ipsec + 0.5) * divphi) * TMath::DegToRad();
489  Double_t cphi = TMath::Cos(phi);
490  Double_t sphi = TMath::Sin(phi);
491  Double_t rproj = point[0] * cphi + point[1] * sphi;
492  if (dir[2] > 0) {
493  ipl = ipln + 1;
494  if (rproj > fRmin[ipln] && rproj < fRmin[ipln + 1]) return 0.0;
495  if (rproj < fRmax[ipln] && rproj > fRmax[ipln + 1]) return 0.0;
496  } else {
497  ipl = ipln - 1;
498  if (rproj < fRmin[ipln] && rproj > fRmin[ipln + 1]) return 0.0;
499  if (rproj > fRmax[ipln] && rproj < fRmax[ipln + 1]) return 0.0;
500  }
501  }
502 
503  Int_t icrossed;
504  icrossed = GetPhiCrossList(point, dir, ipsec, sph, iph, stepmax);
505  Double_t snext;
506  if (TMath::Abs(dir[2]) < TGeoShape::Tolerance()) {
507  if (SliceCrossingInZ(point, dir, icrossed, iph, sph, snext, stepmax)) return snext;
508  if (snext > TGeoShape::Tolerance()) return TGeoShape::Big();
509  return 0.;
510  }
511  if (SliceCrossingIn(point, dir, ipl, icrossed, iph, sph, snext, stepmax)) return snext;
512  if (snext > TGeoShape::Tolerance()) return TGeoShape::Big();
513  return 0.;
514 }
515 
516 ////////////////////////////////////////////////////////////////////////////////
517 /// Locates index IPSEC of the phi sector containing POINT.
518 
519 void TGeoPgon::LocatePhi(const Double_t *point, Int_t &ipsec) const
520 {
521  Double_t phi = TMath::ATan2(point[1], point[0]) * TMath::RadToDeg();
522  while (phi < fPhi1) phi += 360.;
523  ipsec = Int_t(fNedges * (phi - fPhi1) / fDphi); // [0, fNedges-1]
524  if (ipsec > fNedges - 1) ipsec = -1; // in gap
525 }
526 
527 ////////////////////////////////////////////////////////////////////////////////
528 /// Returns lists of PGON phi crossings for a ray starting from POINT.
529 
530 Int_t TGeoPgon::GetPhiCrossList(const Double_t *point, const Double_t *dir, Int_t istart, Double_t *sphi, Int_t *iphi,
531  Double_t stepmax) const
532 {
533  Double_t rxy, phi, cph, sph;
534  Int_t icrossed = 0;
535  if ((1. - TMath::Abs(dir[2])) < 1E-8) {
536  // ray is going parallel with Z
537  iphi[0] = istart;
538  sphi[0] = stepmax;
539  return 1;
540  }
541  Bool_t shootorig = (TMath::Abs(point[0] * dir[1] - point[1] * dir[0]) < 1E-8) ? kTRUE : kFALSE;
542  Double_t divphi = fDphi / fNedges;
543  if (shootorig) {
544  Double_t rdotn = point[0] * dir[0] + point[1] * dir[1];
545  if (rdotn > 0) {
546  sphi[0] = stepmax;
547  iphi[0] = istart;
548  return 1;
549  }
550  sphi[0] = TMath::Sqrt((point[0] * point[0] + point[1] * point[1]) / (1. - dir[2] * dir[2]));
551  iphi[0] = istart;
552  if (sphi[0] > stepmax) {
553  sphi[0] = stepmax;
554  return 1;
555  }
556  phi = TMath::ATan2(dir[1], dir[0]) * TMath::RadToDeg();
557  while (phi < fPhi1) phi += 360.;
558  istart = Int_t((phi - fPhi1) / divphi);
559  if (istart > fNedges - 1) istart = -1;
560  iphi[1] = istart;
561  sphi[1] = stepmax;
562  return 2;
563  }
564  Int_t incsec = Int_t(TMath::Sign(1., point[0] * dir[1] - point[1] * dir[0]));
565  Int_t ist;
566  if (istart < 0)
567  ist = (incsec > 0) ? 0 : fNedges;
568  else
569  ist = (incsec > 0) ? (istart + 1) : istart;
570  Bool_t crossing = kTRUE;
571  Bool_t gapdone = kFALSE;
572  divphi *= TMath::DegToRad();
573  Double_t phi1 = fPhi1 * TMath::DegToRad();
574  while (crossing) {
575  if (istart < 0) gapdone = kTRUE;
576  phi = phi1 + ist * divphi;
577  cph = TMath::Cos(phi);
578  sph = TMath::Sin(phi);
579  crossing = IsCrossingSemiplane(point, dir, cph, sph, sphi[icrossed], rxy);
580  if (!crossing) sphi[icrossed] = stepmax;
581  iphi[icrossed++] = istart;
582  if (crossing) {
583  if (sphi[icrossed - 1] > stepmax) {
584  sphi[icrossed - 1] = stepmax;
585  return icrossed;
586  }
587  if (istart < 0) {
588  istart = (incsec > 0) ? 0 : (fNedges - 1);
589  } else {
590  istart += incsec;
591  if (istart > fNedges - 1)
592  istart = (fDphi < 360.) ? (-1) : 0;
593  else if (istart < 0 && TGeoShape::IsSameWithinTolerance(fDphi, 360))
594  istart = fNedges - 1;
595  }
596  if (istart < 0) {
597  if (gapdone) return icrossed;
598  ist = (incsec > 0) ? 0 : fNedges;
599  } else {
600  ist = (incsec > 0) ? (istart + 1) : istart;
601  }
602  }
603  }
604  return icrossed;
605 }
606 
607 ////////////////////////////////////////////////////////////////////////////////
608 /// Performs ray propagation between Z segments.
609 
610 Bool_t TGeoPgon::SliceCrossingInZ(const Double_t *point, const Double_t *dir, Int_t nphi, Int_t *iphi,
611  Double_t *stepphi, Double_t &snext, Double_t stepmax) const
612 {
613  snext = 0.;
614  if (!nphi) return kFALSE;
615  Int_t i;
616  Double_t rmin, rmax;
617  Double_t apg, bpg;
618  Double_t pt[3];
619  if (iphi[0] < 0 && nphi == 1) return kFALSE;
620  // Get current Z segment
621  Int_t ipl = TMath::BinarySearch(fNz, fZ, point[2]);
622  if (ipl < 0 || ipl == fNz - 1) return kFALSE;
623  if (TMath::Abs(point[2] - fZ[ipl]) < TGeoShape::Tolerance()) {
624  if (ipl < fNz - 2 && TGeoShape::IsSameWithinTolerance(fZ[ipl], fZ[ipl + 1])) {
625  rmin = TMath::Min(fRmin[ipl], fRmin[ipl + 1]);
626  rmax = TMath::Max(fRmax[ipl], fRmax[ipl + 1]);
627  } else if (ipl > 1 && TGeoShape::IsSameWithinTolerance(fZ[ipl], fZ[ipl - 1])) {
628  rmin = TMath::Min(fRmin[ipl], fRmin[ipl + 1]);
629  rmax = TMath::Max(fRmax[ipl], fRmax[ipl + 1]);
630  } else {
631  rmin = fRmin[ipl];
632  rmax = fRmax[ipl];
633  }
634  } else {
635  rmin = Rpg(point[2], ipl, kTRUE, apg, bpg);
636  rmax = Rpg(point[2], ipl, kFALSE, apg, bpg);
637  }
638  Int_t iphcrt;
639  Double_t divphi = TMath::DegToRad() * fDphi / fNedges;
640  Double_t rproj, ndot, dist;
641  Double_t phi1 = fPhi1 * TMath::DegToRad();
642  Double_t cosph, sinph;
643  Double_t snextphi = 0.;
644  Double_t step = 0;
645  Double_t phi;
646  memcpy(pt, point, 3 * sizeof(Double_t));
647  for (iphcrt = 0; iphcrt < nphi; iphcrt++) {
648  if (step > stepmax) {
649  snext = step;
650  return kFALSE;
651  }
652  if (iphi[iphcrt] < 0) {
653  snext = step;
654  return kTRUE;
655  }
656  // check crossing
657  snextphi = stepphi[iphcrt];
658  phi = phi1 + (iphi[iphcrt] + 0.5) * divphi;
659  cosph = TMath::Cos(phi);
660  sinph = TMath::Sin(phi);
661  rproj = pt[0] * cosph + pt[1] * sinph;
662  dist = TGeoShape::Big();
663  ndot = dir[0] * cosph + dir[1] * sinph;
664  if (!TGeoShape::IsSameWithinTolerance(ndot, 0)) {
665  dist = (ndot > 0) ? ((rmax - rproj) / ndot) : ((rmin - rproj) / ndot);
666  if (dist < 0) dist = 0.;
667  }
668  if (dist < (snextphi - step)) {
669  snext = step + dist;
670  if (snext < stepmax) return kTRUE;
671  return kFALSE;
672  }
673  step = snextphi;
674  for (i = 0; i < 3; i++) pt[i] = point[i] + step * dir[i];
675  }
676  snext = step;
677  return kFALSE;
678 }
679 
680 ////////////////////////////////////////////////////////////////////////////////
681 /// Performs ray propagation between Z segments.
682 
683 Bool_t TGeoPgon::SliceCrossingZ(const Double_t *point, const Double_t *dir, Int_t nphi, Int_t *iphi, Double_t *stepphi,
684  Double_t &snext, Double_t stepmax) const
685 {
686  if (!nphi) return kFALSE;
687  Int_t i;
688  Double_t rmin, rmax;
689  Double_t apg, bpg;
690  Double_t pt[3];
691  if (iphi[0] < 0 && nphi == 1) return kFALSE;
692  // Get current Z segment
693  Int_t ipl = TMath::BinarySearch(fNz, fZ, point[2]);
694  if (ipl < 0 || ipl == fNz - 1) return kFALSE;
695  if (TMath::Abs(point[2] - fZ[ipl]) < TGeoShape::Tolerance()) {
696  if (ipl < fNz - 2 && TGeoShape::IsSameWithinTolerance(fZ[ipl], fZ[ipl + 1])) {
697  rmin = TMath::Min(fRmin[ipl], fRmin[ipl + 1]);
698  rmax = TMath::Max(fRmax[ipl], fRmax[ipl + 1]);
699  } else if (ipl > 1 && TGeoShape::IsSameWithinTolerance(fZ[ipl], fZ[ipl - 1])) {
700  rmin = TMath::Min(fRmin[ipl], fRmin[ipl + 1]);
701  rmax = TMath::Max(fRmax[ipl], fRmax[ipl + 1]);
702  } else {
703  rmin = fRmin[ipl];
704  rmax = fRmax[ipl];
705  }
706  } else {
707  rmin = Rpg(point[2], ipl, kTRUE, apg, bpg);
708  rmax = Rpg(point[2], ipl, kFALSE, apg, bpg);
709  }
710  Int_t iphcrt;
711  Double_t divphi = TMath::DegToRad() * fDphi / fNedges;
712  Double_t rproj, ndot, dist;
713  Double_t phi1 = fPhi1 * TMath::DegToRad();
714  Double_t cosph, sinph;
715  Double_t snextphi = 0.;
716  Double_t step = 0;
717  Double_t phi;
718  memcpy(pt, point, 3 * sizeof(Double_t));
719  for (iphcrt = 0; iphcrt < nphi; iphcrt++) {
720  if (step > stepmax) return kFALSE;
721  snextphi = stepphi[iphcrt];
722  if (iphi[iphcrt] < 0) {
723  if (iphcrt == nphi - 1) return kFALSE;
724  if (snextphi > stepmax) return kFALSE;
725  for (i = 0; i < 3; i++) pt[i] = point[i] + snextphi * dir[i];
726  phi = phi1 + (iphi[iphcrt + 1] + 0.5) * divphi;
727  cosph = TMath::Cos(phi);
728  sinph = TMath::Sin(phi);
729  rproj = pt[0] * cosph + pt[1] * sinph;
730  if (rproj < rmin || rproj > rmax) {
731  step = snextphi;
732  continue;
733  }
734  snext = snextphi;
735  return kTRUE;
736  }
737  // check crossing
738  phi = phi1 + (iphi[iphcrt] + 0.5) * divphi;
739  cosph = TMath::Cos(phi);
740  sinph = TMath::Sin(phi);
741  rproj = pt[0] * cosph + pt[1] * sinph;
742  dist = TGeoShape::Big();
743  ndot = dir[0] * cosph + dir[1] * sinph;
744  if (rproj < rmin) {
745  dist = (ndot > 0) ? ((rmin - rproj) / ndot) : TGeoShape::Big();
746  } else {
747  dist = (ndot < 0) ? ((rmax - rproj) / ndot) : TGeoShape::Big();
748  }
749  if (dist < 1E10) {
750  snext = step + dist;
751  if (snext < stepmax) return kTRUE;
752  }
753  step = snextphi;
754  for (i = 0; i < 3; i++) pt[i] = point[i] + step * dir[i];
755  }
756  return kFALSE;
757 }
758 
759 ////////////////////////////////////////////////////////////////////////////////
760 /// Check boundary crossing inside phi slices. Return distance snext to first crossing
761 /// if smaller than stepmax.
762 /// Protection in case point is in phi gap or close to phi boundaries and exiting
763 
764 Bool_t TGeoPgon::SliceCrossingIn(const Double_t *point, const Double_t *dir, Int_t ipl, Int_t nphi, Int_t *iphi,
765  Double_t *stepphi, Double_t &snext, Double_t stepmax) const
766 {
767  snext = 0.;
768  if (!nphi) return kFALSE;
769  Int_t i;
770  Int_t iphstart = 0;
771  Double_t pt[3];
772  if (iphi[0] < 0) {
773  if (stepphi[0] > TGeoShape::Tolerance()) return kFALSE;
774  iphstart = 1;
775  }
776  if (nphi > 1 && iphi[1] < 0 && stepphi[0] < TGeoShape::Tolerance()) {
777  snext = stepphi[0];
778  return kTRUE;
779  }
780  // Get current Z segment
781  Double_t snextphi = 0.;
782  Double_t step = 0;
783  Int_t incseg = (dir[2] > 0) ? 1 : -1; // dir[2] is never 0 here
784  // Compute the projected radius from starting point
785  Int_t iplstart = ipl;
786  Int_t iphcrt = 0;
787  Double_t apr = TGeoShape::Big(), bpr = 0, db = 0;
788  Double_t rpg = 0, rnew = 0, znew = 0;
789  Double_t rpgin = 0, rpgout = 0, apgin = 0, apgout = 0, bpgin = 0, bpgout = 0;
790  Double_t divphi = TMath::DegToRad() * fDphi / fNedges;
791  Double_t phi1 = fPhi1 * TMath::DegToRad();
792  Double_t phi = 0, dz = 0;
793  Double_t cosph = 0, sinph = 0;
794  Double_t distz = 0, distr = 0, din = 0, dout = 0;
795  Double_t invdir = 1. / dir[2];
796  memcpy(pt, point, 3 * sizeof(Double_t));
797  for (iphcrt = iphstart; iphcrt < nphi; iphcrt++) {
798  // check if step to current checked slice is too big
799  if (step > stepmax) {
800  snext = step;
801  return kFALSE;
802  }
803  if (iphi[iphcrt] < 0) {
804  snext = snextphi;
805  return kTRUE;
806  }
807  snextphi = stepphi[iphcrt];
808  phi = phi1 + (iphi[iphcrt] + 0.5) * divphi;
809  cosph = TMath::Cos(phi);
810  sinph = TMath::Sin(phi);
811  Double_t rproj = Rproj(pt[2], pt, dir, cosph, sinph, apr, bpr);
812  // compute distance to next Z plane
813  while (ipl >= 0 && ipl < fNz - 1) {
814  din = dout = TGeoShape::Big();
815  // dist to last boundary of current segment according dir
816  distz = (fZ[ipl + ((1 + incseg) >> 1)] - pt[2]) * invdir;
817  // length of current segment
818  dz = fZ[ipl + 1] - fZ[ipl];
819  if (dz < TGeoShape::Tolerance()) {
820  rnew = apr + bpr * fZ[ipl];
821  rpg = (rnew - fRmin[ipl]) * (rnew - fRmin[ipl + 1]);
822  if (rpg <= 0) din = distz;
823  rpg = (rnew - fRmax[ipl]) * (rnew - fRmax[ipl + 1]);
824  if (rpg <= 0) dout = distz;
825  distr = TMath::Min(din, dout);
826  } else {
827  rpgin = Rpg(pt[2], ipl, kTRUE, apgin, bpgin);
828  db = bpgin - bpr;
829  if (TMath::Abs(db) > TGeoShape::Tolerance()) {
830  znew = (apr - apgin) / db;
831  din = (znew - pt[2]) * invdir;
832  }
833  rpgout = Rpg(pt[2], ipl, kFALSE, apgout, bpgout);
834  db = bpgout - bpr;
835  if (TMath::Abs(db) > TGeoShape::Tolerance()) {
836  znew = (apr - apgout) / db;
837  dout = (znew - pt[2]) * invdir;
838  }
839  // protection for the first segment
840  Double_t dinp = (din > TMath::Abs(snext - TGeoShape::Tolerance())) ? din : TGeoShape::Big();
841  Double_t doutp = (dout > TMath::Abs(snext - TGeoShape::Tolerance())) ? dout : TGeoShape::Big();
842  distr = TMath::Min(dinp, doutp);
843  if (iphcrt == iphstart && ipl == iplstart) {
844  if (rproj < rpgin + 1.E-8) {
845  Double_t ndotd = dir[0] * cosph + dir[1] * sinph + dir[2] * (fRmin[ipl] - fRmin[ipl + 1]) / dz;
846  if (ndotd < 0) {
847  snext = (din < 0) ? step : (step + din);
848  return kTRUE;
849  } else {
850  // Ignore din
851  din = -TGeoShape::Big();
852  }
853  distr = TMath::Max(din, dout);
854  if (distr < TGeoShape::Tolerance()) distr = TGeoShape::Big();
855  } else if (rproj > rpgout - 1.E-8) {
856  Double_t ndotd = dir[0] * cosph + dir[1] * sinph + dir[2] * (fRmax[ipl] - fRmax[ipl + 1]) / dz;
857  if (ndotd > 0) {
858  snext = (dout < 0) ? step : (step + dout);
859  return kTRUE;
860  } else {
861  // Ignore dout
862  dout = -TGeoShape::Big();
863  }
864  distr = TMath::Max(din, dout);
865  if (distr < TGeoShape::Tolerance()) distr = TGeoShape::Big();
866  }
867  }
868  }
869  if (distr < snext - TGeoShape::Tolerance()) distr = TGeoShape::Big();
870  if (snextphi < step + TMath::Min(distz, distr)) {
871  for (i = 0; i < 3; i++) pt[i] = point[i] + snextphi * dir[i];
872  step = snextphi;
873  snext = 0.0;
874  break;
875  }
876  if (distr <= distz + TGeoShape::Tolerance()) {
877  step += distr;
878  snext = step;
879  return (step > stepmax) ? kFALSE : kTRUE;
880  }
881  // we have crossed a Z boundary
882  snext = distz;
883  if ((ipl + incseg < 0) || (ipl + incseg > fNz - 2)) {
884  // it was the last boundary
885  step += distz;
886  snext = step;
887  return (step > stepmax) ? kFALSE : kTRUE;
888  }
889  ipl += incseg;
890  } // end loop Z
891  } // end loop phi
892  snext = TGeoShape::Big();
893  return kFALSE;
894 }
895 
896 ////////////////////////////////////////////////////////////////////////////////
897 /// Check boundary crossing inside phi slices. Return distance snext to first crossing
898 /// if smaller than stepmax.
899 
900 Bool_t TGeoPgon::SliceCrossing(const Double_t *point, const Double_t *dir, Int_t nphi, Int_t *iphi, Double_t *stepphi,
901  Double_t &snext, Double_t stepmax) const
902 {
903  if (!nphi) return kFALSE;
904  Int_t i;
905  Double_t pt[3];
906  if (iphi[0] < 0 && nphi == 1) return kFALSE;
907 
908  Double_t snextphi = 0.;
909  Double_t step = 0;
910  // Get current Z segment
911  Int_t incseg = (dir[2] > 0) ? 1 : -1; // dir[2] is never 0 here
912  Int_t ipl = TMath::BinarySearch(fNz, fZ, point[2]);
913  if (ipl < 0) {
914  ipl = 0; // this should never happen
915  if (incseg < 0) return kFALSE;
916  } else {
917  if (ipl == fNz - 1) {
918  ipl = fNz - 2; // nor this
919  if (incseg > 0) return kFALSE;
920  } else {
921  if (TMath::Abs(point[2] - fZ[ipl]) < TGeoShape::Tolerance()) {
922  // we are at the sector edge, but never inside the pgon
923  if ((ipl + incseg) < 0 || (ipl + incseg) > fNz - 1) return kFALSE;
924  if (TGeoShape::IsSameWithinTolerance(fZ[ipl], fZ[ipl + incseg])) ipl += incseg;
925  // move to next clean segment if downwards
926  if (incseg < 0) {
927  if (TGeoShape::IsSameWithinTolerance(fZ[ipl], fZ[ipl + 1])) ipl--;
928  }
929  }
930  }
931  }
932  // Compute the projected radius from starting point
933  Int_t iphcrt;
934  Double_t apg, bpg;
935  Double_t rpgin;
936  Double_t rpgout;
937  Double_t divphi = TMath::DegToRad() * fDphi / fNedges;
938  Double_t phi1 = fPhi1 * TMath::DegToRad();
939  Double_t phi;
940  Double_t cosph, sinph;
941  Double_t rproj;
942  memcpy(pt, point, 3 * sizeof(Double_t));
943  for (iphcrt = 0; iphcrt < nphi; iphcrt++) {
944  // check if step to current checked slice is too big
945  if (step > stepmax) return kFALSE;
946  // jump over the dead sector
947  snextphi = stepphi[iphcrt];
948  if (iphi[iphcrt] < 0) {
949  if (iphcrt == nphi - 1) return kFALSE;
950  if (snextphi > stepmax) return kFALSE;
951  for (i = 0; i < 3; i++) pt[i] = point[i] + snextphi * dir[i];
952  // we have a new z, so check again iz
953  if (incseg > 0) {
954  // loop z planes
955  while (pt[2] > fZ[ipl + 1]) {
956  ipl++;
957  if (ipl > fNz - 2) return kFALSE;
958  }
959  } else {
960  while (pt[2] < fZ[ipl]) {
961  ipl--;
962  if (ipl < 0) return kFALSE;
963  }
964  }
965  // check if we have a crossing when entering new sector
966  rpgin = Rpg(pt[2], ipl, kTRUE, apg, bpg);
967  rpgout = Rpg(pt[2], ipl, kFALSE, apg, bpg);
968  phi = phi1 + (iphi[iphcrt + 1] + 0.5) * divphi;
969  cosph = TMath::Cos(phi);
970  sinph = TMath::Sin(phi);
971 
972  rproj = pt[0] * cosph + pt[1] * sinph;
973  if (rproj < rpgin || rproj > rpgout) {
974  step = snextphi;
975  continue;
976  }
977  snext = snextphi;
978  return kTRUE;
979  }
980  if (IsCrossingSlice(point, dir, iphi[iphcrt], step, ipl, snext, TMath::Min(snextphi, stepmax))) return kTRUE;
981  step = snextphi;
982  }
983  return kFALSE;
984 }
985 
986 ////////////////////////////////////////////////////////////////////////////////
987 /// Check crossing of a given pgon slice, from a starting point inside the slice
988 
989 Bool_t TGeoPgon::IsCrossingSlice(const Double_t *point, const Double_t *dir, Int_t iphi, Double_t sstart, Int_t &ipl,
990  Double_t &snext, Double_t stepmax) const
991 {
992  if (ipl < 0 || ipl > fNz - 2) return kFALSE;
993  if (sstart > stepmax) return kFALSE;
994  Double_t pt[3];
995  memcpy(pt, point, 3 * sizeof(Double_t));
996  if (sstart > 0)
997  for (Int_t i = 0; i < 3; i++) pt[i] += sstart * dir[i];
998  stepmax -= sstart;
999  Double_t step;
1000  Int_t incseg = (dir[2] > 0) ? 1 : -1;
1001  Double_t invdir = 1. / dir[2];
1002  Double_t divphi = TMath::DegToRad() * fDphi / fNedges;
1003  Double_t phi = fPhi1 * TMath::DegToRad() + (iphi + 0.5) * divphi;
1004  Double_t cphi = TMath::Cos(phi);
1005  Double_t sphi = TMath::Sin(phi);
1006  Double_t apr = TGeoShape::Big();
1007  Double_t bpr = 0.;
1008  Rproj(pt[2], point, dir, cphi, sphi, apr, bpr);
1009  Double_t dz;
1010  // loop segments
1011  Int_t icrtseg = ipl;
1012  Int_t isegstart = ipl;
1013  Int_t iseglast = (incseg > 0) ? (fNz - 1) : -1;
1014  Double_t din, dout, rdot, rnew, rpg, apg, bpg, db, znew;
1015 
1016  for (ipl = isegstart; ipl != iseglast; ipl += incseg) {
1017  step = (fZ[ipl + 1 - ((1 + incseg) >> 1)] - pt[2]) * invdir;
1018  if (step > 0) {
1019  if (step > stepmax) {
1020  ipl = icrtseg;
1021  return kFALSE;
1022  }
1023  icrtseg = ipl;
1024  }
1025  din = dout = TGeoShape::Big();
1026  dz = fZ[ipl + 1] - fZ[ipl];
1027 
1028  // rdot = (rproj-fRmin[ipl])*dz - (pt[2]-fZ[ipl])*(fRmin[ipl+1]-fRmin[ipl]);
1029  if (TGeoShape::IsSameWithinTolerance(dz, 0))
1030  rdot = dir[2] * TMath::Sign(1., fRmin[ipl] - fRmin[ipl + 1]);
1031  else
1032  rdot = dir[0] * cphi + dir[1] * sphi + dir[2] * (fRmin[ipl] - fRmin[ipl + 1]) / dz;
1033  if (rdot > 0) {
1034  // inner surface visible ->check crossing
1035  // printf(" inner visible\n");
1036  if (TGeoShape::IsSameWithinTolerance(dz, 0)) {
1037  rnew = apr + bpr * fZ[ipl];
1038  rpg = (rnew - fRmin[ipl]) * (rnew - fRmin[ipl + 1]);
1039  if (rpg <= 0) din = (fZ[ipl] - pt[2]) * invdir;
1040  } else {
1041  rpg = Rpg(pt[2], ipl, kTRUE, apg, bpg);
1042  db = bpg - bpr;
1043  if (!TGeoShape::IsSameWithinTolerance(db, 0)) {
1044  znew = (apr - apg) / db;
1045  if (znew > fZ[ipl] && znew < fZ[ipl + 1]) {
1046  din = (znew - pt[2]) * invdir;
1047  if (din < 0) din = TGeoShape::Big();
1048  }
1049  }
1050  }
1051  }
1052  // printf(" din=%f\n", din);
1053  // rdot = (rproj-fRmax[ipl])*dz - (pt[2]-fZ[ipl])*(fRmax[ipl+1]-fRmax[ipl]);
1054  if (TGeoShape::IsSameWithinTolerance(dz, 0))
1055  rdot = dir[2] * TMath::Sign(1., fRmax[ipl] - fRmax[ipl + 1]);
1056  else
1057  rdot = dir[0] * cphi + dir[1] * sphi + dir[2] * (fRmax[ipl] - fRmax[ipl + 1]) / dz;
1058  if (rdot < 0) {
1059  // printf(" outer visible\n");
1060  // outer surface visible ->check crossing
1061  if (TGeoShape::IsSameWithinTolerance(dz, 0)) {
1062  rnew = apr + bpr * fZ[ipl];
1063  rpg = (rnew - fRmax[ipl]) * (rnew - fRmax[ipl + 1]);
1064  if (rpg <= 0) dout = (fZ[ipl] - pt[2]) * invdir;
1065  } else {
1066  rpg = Rpg(pt[2], ipl, kFALSE, apg, bpg);
1067  db = bpg - bpr;
1068  if (!TGeoShape::IsSameWithinTolerance(db, 0)) {
1069  znew = (apr - apg) / db;
1070  if (znew > fZ[ipl] && znew < fZ[ipl + 1]) dout = (znew - pt[2]) * invdir;
1071  if (dout < 0) dout = TGeoShape::Big();
1072  }
1073  }
1074  }
1075  // printf(" dout=%f\n", dout);
1076  step = TMath::Min(din, dout);
1077  if (step < 1E10) {
1078  // there is a crossing within this segment
1079  if (step > stepmax) {
1080  ipl = icrtseg;
1081  return kFALSE;
1082  }
1083  snext = sstart + step;
1084  return kTRUE;
1085  }
1086  }
1087  ipl = icrtseg;
1088  return kFALSE;
1089 }
1090 
1091 ////////////////////////////////////////////////////////////////////////////////
1092 /// Compute distance from outside point to surface of the polygone
1093 
1094 Double_t TGeoPgon::DistFromOutside(const Double_t *point, const Double_t *dir, Int_t iact, Double_t step,
1095  Double_t *safe) const
1096 {
1097  if (iact < 3 && safe) {
1098  *safe = Safety(point, kFALSE);
1099  if (iact == 0) return TGeoShape::Big(); // just safety computed
1100  if (iact == 1 && step < *safe) return TGeoShape::Big(); // safety mode
1101  }
1102  // Check if the bounding box is crossed within the requested distance
1103  Double_t sdist = TGeoBBox::DistFromOutside(point, dir, fDX, fDY, fDZ, fOrigin, step);
1104  if (sdist >= step) return TGeoShape::Big();
1105  // Protection for points on last Z sections
1106  if (dir[2] <= 0 && TMath::Abs(point[2] - fZ[0]) < TGeoShape::Tolerance()) return TGeoShape::Big();
1107  if (dir[2] >= 0 && TMath::Abs(point[2] - fZ[fNz - 1]) < TGeoShape::Tolerance()) return TGeoShape::Big();
1108  // copy the current point
1109  Double_t pt[3];
1110  memcpy(pt, point, 3 * sizeof(Double_t));
1111  // find current Z section
1112  Int_t ipl;
1113  Int_t i, ipsec;
1114  ipl = TMath::BinarySearch(fNz, fZ, pt[2]);
1115 
1116  Double_t divphi = fDphi / fNedges;
1117  // check if ray may intersect outer cylinder
1118  Double_t snext = 0.;
1119  Double_t stepmax = step;
1120  Double_t rpr, snewcross;
1121  Double_t r2 = pt[0] * pt[0] + pt[1] * pt[1];
1122  Double_t radmax = fRmax[TMath::LocMax(fNz, fRmax)];
1123  radmax = radmax / TMath::Cos(0.5 * divphi * TMath::DegToRad());
1124  radmax += 1E-8;
1125  if (r2 > (radmax * radmax) || pt[2] < fZ[0] || pt[2] > fZ[fNz - 1]) {
1126  pt[2] -= 0.5 * (fZ[0] + fZ[fNz - 1]);
1127  snext = TGeoTube::DistFromOutsideS(pt, dir, 0., radmax, 0.5 * (fZ[fNz - 1] - fZ[0]));
1128  if (snext > 1E10) return TGeoShape::Big();
1129  if (snext > stepmax) return TGeoShape::Big();
1130  stepmax -= snext;
1131  pt[2] = point[2];
1132  for (i = 0; i < 3; i++) pt[i] += snext * dir[i];
1133  Bool_t checkz = (ipl < 0 && TMath::Abs(pt[2] - fZ[0]) < 1E-8) ? kTRUE : kFALSE;
1134  if (!checkz) checkz = (ipl == fNz - 1 && TMath::Abs(pt[2] - fZ[fNz - 1]) < 1E-8) ? kTRUE : kFALSE;
1135  if (checkz) {
1136  Double_t rmin, rmax;
1137  if (ipl < 0) {
1138  rmin = fRmin[0];
1139  rmax = fRmax[0];
1140  } else {
1141  rmin = fRmin[fNz - 1];
1142  rmax = fRmax[fNz - 1];
1143  }
1144  Double_t phi = TMath::ATan2(pt[1], pt[0]) * TMath::RadToDeg();
1145  while (phi < fPhi1) phi += 360.0;
1146  Double_t ddp = phi - fPhi1;
1147  if (ddp <= fDphi) {
1148  ipsec = Int_t(ddp / divphi);
1149  Double_t ph0 = (fPhi1 + divphi * (ipsec + 0.5)) * TMath::DegToRad();
1150  rpr = pt[0] * TMath::Cos(ph0) + pt[1] * TMath::Sin(ph0);
1151  if (rpr >= rmin && rpr <= rmax) return snext;
1152  }
1153  }
1154  }
1155  if (!fThreadSize) ((TGeoPgon *)this)->CreateThreadData(1);
1156  ThreadData_t &td = GetThreadData();
1157  Double_t *sph = td.fDblBuffer;
1158  Int_t *iph = td.fIntBuffer;
1159  Int_t icrossed;
1160  // locate current phi sector [0,fNedges-1]; -1 for dead region
1161  // if ray is perpendicular to Z, solve this particular case
1162  if (TMath::Abs(dir[2]) < TGeoShape::Tolerance()) {
1163  LocatePhi(pt, ipsec);
1164  icrossed = GetPhiCrossList(pt, dir, ipsec, sph, iph, stepmax);
1165  if (SliceCrossingZ(pt, dir, icrossed, iph, sph, snewcross, stepmax)) return (snext + snewcross);
1166  return TGeoShape::Big();
1167  }
1168  // Locate phi and get the phi crossing list
1169  divphi *= TMath::DegToRad();
1170  Bool_t inphi = kTRUE;
1171  Double_t ph = TMath::ATan2(pt[1], pt[0]) * TMath::RadToDeg();
1172  while (ph < fPhi1) ph += 360.;
1173  ipsec = Int_t(fNedges * (ph - fPhi1) / fDphi); // [0, fNedges-1]
1174  if (ipsec > fNedges - 1) ipsec = -1; // in gap
1175  Double_t phim = fPhi1 + 0.5 * fDphi;
1176  Double_t ddp = TMath::Abs(ph - phim);
1177  if (fDphi < 360.0) {
1178  inphi = (ddp < 0.5 * fDphi + TGeoShape::Tolerance()) ? kTRUE : kFALSE;
1179  }
1180  ipl = TMath::BinarySearch(fNz, fZ, pt[2]);
1181  if (ipl < 0) ipl = 0;
1182  if (ipl == fNz - 1) ipl--;
1183  Bool_t inz = kTRUE;
1184  if (pt[2] > fZ[fNz - 1] + TGeoShape::Tolerance()) inz = kFALSE;
1185  if (pt[2] < fZ[0] - TGeoShape::Tolerance()) inz = kFALSE;
1186  Bool_t onphi = kFALSE;
1187  if (inphi && inz) {
1188  Bool_t done = kFALSE;
1189  Double_t dz = fZ[ipl + 1] - fZ[ipl];
1190  Double_t phi = fPhi1 * TMath::DegToRad() + (ipsec + 0.5) * divphi;
1191  Double_t cphi = TMath::Cos(phi);
1192  Double_t sphi = TMath::Sin(phi);
1193  Double_t rproj = pt[0] * cphi + pt[1] * sphi;
1194  if (TGeoShape::IsSameWithinTolerance(dz, 0)) {
1195  if (rproj < fRmin[ipl] && rproj > fRmin[ipl + 1] && dir[2] > 0) return 0.0;
1196  if (rproj > fRmin[ipl] && rproj < fRmin[ipl + 1] && dir[2] < 0) return 0.0;
1197  if (rproj > fRmax[ipl] && rproj < fRmax[ipl + 1] && dir[2] > 0) return 0.0;
1198  if (rproj < fRmax[ipl] && rproj > fRmax[ipl + 1] && dir[2] < 0) return 0.0;
1199  done = kTRUE;
1200  }
1201  if (!done) {
1202  Double_t apgout, bpgout;
1203  Double_t rpgout = Rpg(pt[2], ipl, kFALSE, apgout, bpgout);
1204  if (rproj < rpgout + 1.E-8) {
1205  Double_t apgin, bpgin;
1206  Double_t rpgin = Rpg(pt[2], ipl, kTRUE, apgin, bpgin);
1207  if (rproj > rpgin - 1.E-8) {
1208  Double_t safrmin = rproj - rpgin;
1209  Double_t safrmax = rpgout - rproj;
1210  Double_t safz = TMath::Min(pt[2] - fZ[ipl], fZ[ipl + 1] - pt[2]);
1211  Double_t safphi = TGeoShape::Big();
1212  if (fDphi < 360) {
1213  safphi = rproj * TMath::Sin((ddp - 0.5 * fDphi) * TMath::DegToRad());
1214  safphi = TMath::Abs(safphi);
1215  }
1216  // printf("inside pgon: safrmin=%f, safrmax=%f, safphi=%f,
1217  // safz=%f\n",safrmin,safrmax,safphi,safz);
1218  Double_t dzinv = 1. / dz;
1219  if (safrmin < safz && safrmin < safrmax && safrmin < safphi) {
1220  // on inner boundary
1221  Double_t ndotd = dir[0] * cphi + dir[1] * sphi + dir[2] * (fRmin[ipl] - fRmin[ipl + 1]) * dzinv;
1222  // printf(" - inner ndotd=%f (>0 ->0)\n",ndotd);
1223  if (ndotd > 0) return snext;
1224  done = kTRUE;
1225  }
1226  if (!done && safrmax < safz && safrmax < safphi) {
1227  Double_t ndotd = dir[0] * cphi + dir[1] * sphi + dir[2] * (fRmax[ipl] - fRmax[ipl + 1]) * dzinv;
1228  // printf(" - outer ndotd=%f (<0 ->0)\n",ndotd);
1229  if (ndotd < 0) return snext;
1230  done = kTRUE;
1231  }
1232  if (!done && safz < safphi) {
1233  done = kTRUE;
1234  Int_t iplc = ipl;
1235  if (TMath::Abs(pt[2] - fZ[ipl]) > TMath::Abs(fZ[ipl + 1] - pt[2])) iplc++;
1236  if (iplc == 0 || iplc == fNz - 1) {
1237  if (pt[2] * dir[2] < 0) return snext;
1238  return TGeoShape::Big();
1239  } else {
1240  if (TGeoShape::IsSameWithinTolerance(fZ[iplc], fZ[iplc + 1])) {
1241  if (dir[2] > 0) {
1242  if (rproj < fRmin[iplc] && rproj > fRmin[iplc + 1]) return snext;
1243  if (rproj > fRmax[iplc] && rproj < fRmax[iplc + 1]) return snext;
1244  } else {
1245  if (rproj > fRmin[iplc] && rproj < fRmin[iplc + 1]) return snext;
1246  if (rproj < fRmax[iplc] && rproj > fRmax[iplc + 1]) return snext;
1247  }
1248  } else if (TGeoShape::IsSameWithinTolerance(fZ[iplc], fZ[iplc - 1])) {
1249  if (dir[2] > 0) {
1250  if (rproj < fRmin[iplc - 1] && rproj > fRmin[iplc]) return snext;
1251  if (rproj > fRmax[iplc - 1] && rproj < fRmax[iplc]) return snext;
1252  } else {
1253  if (rproj > fRmin[iplc - 1] && rproj < fRmin[iplc]) return snext;
1254  if (rproj < fRmax[iplc - 1] && rproj > fRmax[iplc]) return snext;
1255  }
1256  }
1257  }
1258  }
1259  if (!done) {
1260  // point on phi boundary
1261  onphi = kTRUE;
1262  }
1263  }
1264  }
1265  }
1266  }
1267  icrossed = GetPhiCrossList(pt, dir, ipsec, sph, iph, stepmax);
1268  if (onphi) {
1269  if (!icrossed) return snext;
1270  if (iph[0] < 0 && sph[0] < TGeoShape::Tolerance()) return (snext + sph[0]);
1271  if (iph[0] >= 0 && sph[0] > 1.E-8) return snext;
1272  }
1273  // Fire-up slice crossing algorithm
1274  if (SliceCrossing(pt, dir, icrossed, iph, sph, snewcross, stepmax)) {
1275  snext += snewcross;
1276  return snext;
1277  }
1278  return TGeoShape::Big();
1279 }
1280 
1281 ////////////////////////////////////////////////////////////////////////////////
1282 /// compute closest distance from point px,py to each corner
1283 
1284 Int_t TGeoPgon::DistancetoPrimitive(Int_t px, Int_t py)
1285 {
1286  Int_t n = fNedges + 1;
1287  const Int_t numPoints = 2 * n * fNz;
1288  return ShapeDistancetoPrimitive(numPoints, px, py);
1289 }
1290 
1291 ////////////////////////////////////////////////////////////////////////////////
1292 /// Divide this polygone shape belonging to volume "voldiv" into ndiv volumes
1293 /// called divname, from start position with the given step. Returns pointer
1294 /// to created division cell volume in case of Z divisions. Phi divisions are
1295 /// allowed only if nedges%ndiv=0 and create polygone "segments" with nedges/ndiv edges.
1296 /// Z divisions can be performed if the divided range is in between two consecutive Z planes.
1297 /// In case a wrong division axis is supplied, returns pointer to volume that was divided.
1298 
1299 TGeoVolume *TGeoPgon::Divide(TGeoVolume *voldiv, const char *divname, Int_t iaxis, Int_t ndiv, Double_t start,
1300  Double_t step)
1301 {
1302  // printf("Dividing %s : nz=%d nedges=%d phi1=%g dphi=%g (ndiv=%d iaxis=%d start=%g step=%g)\n",
1303  // voldiv->GetName(), fNz, fNedges, fPhi1, fDphi, ndiv, iaxis, start, step);
1304  TGeoShape *shape; //--- shape to be created
1305  TGeoVolume *vol; //--- division volume to be created
1306  TGeoVolumeMulti *vmulti; //--- generic divided volume
1307  TGeoPatternFinder *finder; //--- finder to be attached
1308  TString opt = ""; //--- option to be attached
1309  Int_t nedges = fNedges;
1310  Double_t zmin = start;
1311  Double_t zmax = start + ndiv * step;
1312  Int_t isect = -1;
1313  Int_t is, id, ipl;
1314  switch (iaxis) {
1315  case 1: //--- R division
1316  Error("Divide", "makes no sense dividing a pgon on radius");
1317  return 0;
1318  case 2: //--- Phi division
1319  if (fNedges % ndiv) {
1320  Error("Divide", "ndiv should divide number of pgon edges");
1321  return 0;
1322  }
1323  nedges = fNedges / ndiv;
1324  finder = new TGeoPatternCylPhi(voldiv, ndiv, start, start + ndiv * step);
1325  vmulti = gGeoManager->MakeVolumeMulti(divname, voldiv->GetMedium());
1326  voldiv->SetFinder(finder);
1327  finder->SetDivIndex(voldiv->GetNdaughters());
1328  shape = new TGeoPgon(-step / 2, step, nedges, fNz);
1329  vol = new TGeoVolume(divname, shape, voldiv->GetMedium());
1330  vmulti->AddVolume(vol);
1331  for (is = 0; is < fNz; is++) ((TGeoPgon *)shape)->DefineSection(is, fZ[is], fRmin[is], fRmax[is]);
1332  opt = "Phi";
1333  for (id = 0; id < ndiv; id++) {
1334  voldiv->AddNodeOffset(vol, id, start + id * step + step / 2, opt.Data());
1335  ((TGeoNodeOffset *)voldiv->GetNodes()->At(voldiv->GetNdaughters() - 1))->SetFinder(finder);
1336  }
1337  return vmulti;
1338  case 3: // --- Z division
1339  // find start plane
1340  for (ipl = 0; ipl < fNz - 1; ipl++) {
1341  if (start < fZ[ipl])
1342  continue;
1343  else {
1344  if ((start + ndiv * step) > fZ[ipl + 1]) continue;
1345  }
1346  isect = ipl;
1347  zmin = fZ[isect];
1348  zmax = fZ[isect + 1];
1349  break;
1350  }
1351  if (isect < 0) {
1352  Error("Divide", "cannot divide pcon on Z if divided region is not between 2 consecutive planes");
1353  return 0;
1354  }
1355  finder = new TGeoPatternZ(voldiv, ndiv, start, start + ndiv * step);
1356  vmulti = gGeoManager->MakeVolumeMulti(divname, voldiv->GetMedium());
1357  voldiv->SetFinder(finder);
1358  finder->SetDivIndex(voldiv->GetNdaughters());
1359  opt = "Z";
1360  for (id = 0; id < ndiv; id++) {
1361  Double_t z1 = start + id * step;
1362  Double_t z2 = start + (id + 1) * step;
1363  Double_t rmin1 = (fRmin[isect] * (zmax - z1) - fRmin[isect + 1] * (zmin - z1)) / (zmax - zmin);
1364  Double_t rmax1 = (fRmax[isect] * (zmax - z1) - fRmax[isect + 1] * (zmin - z1)) / (zmax - zmin);
1365  Double_t rmin2 = (fRmin[isect] * (zmax - z2) - fRmin[isect + 1] * (zmin - z2)) / (zmax - zmin);
1366  Double_t rmax2 = (fRmax[isect] * (zmax - z2) - fRmax[isect + 1] * (zmin - z2)) / (zmax - zmin);
1367  shape = new TGeoPgon(fPhi1, fDphi, nedges, 2);
1368  ((TGeoPgon *)shape)->DefineSection(0, -step / 2, rmin1, rmax1);
1369  ((TGeoPgon *)shape)->DefineSection(1, step / 2, rmin2, rmax2);
1370  vol = new TGeoVolume(divname, shape, voldiv->GetMedium());
1371  vmulti->AddVolume(vol);
1372  voldiv->AddNodeOffset(vol, id, start + id * step + step / 2, opt.Data());
1373  ((TGeoNodeOffset *)voldiv->GetNodes()->At(voldiv->GetNdaughters() - 1))->SetFinder(finder);
1374  }
1375  return vmulti;
1376  default: Error("Divide", "Wrong axis type for division"); return 0;
1377  }
1378 }
1379 
1380 ////////////////////////////////////////////////////////////////////////////////
1381 /// Fill vector param[4] with the bounding cylinder parameters. The order
1382 /// is the following : Rmin, Rmax, Phi1, Phi2
1383 
1384 void TGeoPgon::GetBoundingCylinder(Double_t *param) const
1385 {
1386  param[0] = fRmin[0]; // Rmin
1387  param[1] = fRmax[0]; // Rmax
1388  for (Int_t i = 1; i < fNz; i++) {
1389  if (fRmin[i] < param[0]) param[0] = fRmin[i];
1390  if (fRmax[i] > param[1]) param[1] = fRmax[i];
1391  }
1392  Double_t divphi = fDphi / fNedges;
1393  param[1] /= TMath::Cos(0.5 * divphi * TMath::DegToRad());
1394  param[0] *= param[0];
1395  param[1] *= param[1];
1396  if (TGeoShape::IsSameWithinTolerance(fDphi, 360)) {
1397  param[2] = 0.;
1398  param[3] = 360.;
1399  return;
1400  }
1401  param[2] = (fPhi1 < 0) ? (fPhi1 + 360.) : fPhi1; // Phi1
1402  param[3] = param[2] + fDphi; // Phi2
1403 }
1404 
1405 ////////////////////////////////////////////////////////////////////////////////
1406 /// Inspect the PGON parameters.
1407 
1408 void TGeoPgon::InspectShape() const
1409 {
1410  printf("*** Shape %s: TGeoPgon ***\n", GetName());
1411  printf(" Nedges = %i\n", fNedges);
1412  TGeoPcon::InspectShape();
1413 }
1414 
1415 ////////////////////////////////////////////////////////////////////////////////
1416 /// Creates a TBuffer3D describing *this* shape.
1417 /// Coordinates are in local reference frame.
1418 
1419 TBuffer3D *TGeoPgon::MakeBuffer3D() const
1420 {
1421  const Int_t n = GetNsegments() + 1;
1422  Int_t nz = GetNz();
1423  if (nz < 2) return 0;
1424  Int_t nbPnts = nz * 2 * n;
1425  if (nbPnts <= 0) return 0;
1426  Double_t dphi = GetDphi();
1427  Bool_t specialCase = kFALSE;
1428  if (TGeoShape::IsSameWithinTolerance(dphi, 360)) specialCase = kTRUE;
1429  Int_t nbSegs = 4 * (nz * n - 1 + (specialCase == kTRUE));
1430  Int_t nbPols = 2 * (nz * n - 1 + (specialCase == kTRUE));
1431 
1432  TBuffer3D *buff =
1433  new TBuffer3D(TBuffer3DTypes::kGeneric, nbPnts, 3 * nbPnts, nbSegs, 3 * nbSegs, nbPols, 6 * nbPols);
1434  if (buff) {
1435  SetPoints(buff->fPnts);
1436  SetSegsAndPols(*buff);
1437  }
1438 
1439  return buff;
1440 }
1441 
1442 ////////////////////////////////////////////////////////////////////////////////
1443 /// Fill TBuffer3D structure for segments and polygons.
1444 
1445 void TGeoPgon::SetSegsAndPols(TBuffer3D &buff) const
1446 {
1447  Int_t i, j;
1448  const Int_t n = GetNsegments() + 1;
1449  Int_t nz = GetNz();
1450  if (nz < 2) return;
1451  Int_t nbPnts = nz * 2 * n;
1452  if (nbPnts <= 0) return;
1453  Double_t dphi = GetDphi();
1454  Bool_t specialCase = kFALSE;
1455  if (TGeoShape::IsSameWithinTolerance(dphi, 360)) specialCase = kTRUE;
1456  Int_t c = GetBasicColor();
1457 
1458  Int_t indx, indx2, k;
1459  indx = indx2 = 0;
1460 
1461  // inside & outside circles, number of segments: 2*nz*(n-1)
1462  // special case number of segments: 2*nz*n
1463  for (i = 0; i < nz * 2; i++) {
1464  indx2 = i * n;
1465  for (j = 1; j < n; j++) {
1466  buff.fSegs[indx++] = c;
1467  buff.fSegs[indx++] = indx2 + j - 1;
1468  buff.fSegs[indx++] = indx2 + j;
1469  }
1470  if (specialCase) {
1471  buff.fSegs[indx++] = c;
1472  buff.fSegs[indx++] = indx2 + j - 1;
1473  buff.fSegs[indx++] = indx2;
1474  }
1475  }
1476 
1477  // bottom & top lines, number of segments: 2*n
1478  for (i = 0; i < 2; i++) {
1479  indx2 = i * (nz - 1) * 2 * n;
1480  for (j = 0; j < n; j++) {
1481  buff.fSegs[indx++] = c;
1482  buff.fSegs[indx++] = indx2 + j;
1483  buff.fSegs[indx++] = indx2 + n + j;
1484  }
1485  }
1486 
1487  // inside & outside cylinders, number of segments: 2*(nz-1)*n
1488  for (i = 0; i < (nz - 1); i++) {
1489  // inside cylinder
1490  indx2 = i * n * 2;
1491  for (j = 0; j < n; j++) {
1492  buff.fSegs[indx++] = c + 2;
1493  buff.fSegs[indx++] = indx2 + j;
1494  buff.fSegs[indx++] = indx2 + n * 2 + j;
1495  }
1496  // outside cylinder
1497  indx2 = i * n * 2 + n;
1498  for (j = 0; j < n; j++) {
1499  buff.fSegs[indx++] = c + 3;
1500  buff.fSegs[indx++] = indx2 + j;
1501  buff.fSegs[indx++] = indx2 + n * 2 + j;
1502  }
1503  }
1504 
1505  // left & right sections, number of segments: 2*(nz-2)
1506  // special case number of segments: 0
1507  if (!specialCase) {
1508  for (i = 1; i < (nz - 1); i++) {
1509  for (j = 0; j < 2; j++) {
1510  buff.fSegs[indx++] = c;
1511  buff.fSegs[indx++] = 2 * i * n + j * (n - 1);
1512  buff.fSegs[indx++] = (2 * i + 1) * n + j * (n - 1);
1513  }
1514  }
1515  }
1516 
1517  Int_t m = n - 1 + (specialCase == kTRUE);
1518  indx = 0;
1519 
1520  // bottom & top, number of polygons: 2*(n-1)
1521  // special case number of polygons: 2*n
1522  i = 0;
1523  for (j = 0; j < n - 1; j++) {
1524  buff.fPols[indx++] = c + 3;
1525  buff.fPols[indx++] = 4;
1526  buff.fPols[indx++] = 2 * nz * m + i * n + j;
1527  buff.fPols[indx++] = i * (nz * 2 - 2) * m + m + j;
1528  buff.fPols[indx++] = 2 * nz * m + i * n + j + 1;
1529  buff.fPols[indx++] = i * (nz * 2 - 2) * m + j;
1530  }
1531  if (specialCase) {
1532  buff.fPols[indx++] = c + 3;
1533  buff.fPols[indx++] = 4;
1534  buff.fPols[indx++] = 2 * nz * m + i * n + j;
1535  buff.fPols[indx++] = i * (nz * 2 - 2) * m + m + j;
1536  buff.fPols[indx++] = 2 * nz * m + i * n;
1537  buff.fPols[indx++] = i * (nz * 2 - 2) * m + j;
1538  }
1539  i = 1;
1540  for (j = 0; j < n - 1; j++) {
1541  buff.fPols[indx++] = c + 3;
1542  buff.fPols[indx++] = 4;
1543  buff.fPols[indx++] = i * (nz * 2 - 2) * m + j;
1544  buff.fPols[indx++] = 2 * nz * m + i * n + j + 1;
1545  buff.fPols[indx++] = i * (nz * 2 - 2) * m + m + j;
1546  buff.fPols[indx++] = 2 * nz * m + i * n + j;
1547  }
1548  if (specialCase) {
1549  buff.fPols[indx++] = c + 3;
1550  buff.fPols[indx++] = 4;
1551  buff.fPols[indx++] = i * (nz * 2 - 2) * m + j;
1552  buff.fPols[indx++] = 2 * nz * m + i * n;
1553  buff.fPols[indx++] = i * (nz * 2 - 2) * m + m + j;
1554  buff.fPols[indx++] = 2 * nz * m + i * n + j;
1555  }
1556 
1557  // inside & outside, number of polygons: (nz-1)*2*(n-1)
1558  for (k = 0; k < (nz - 1); k++) {
1559  i = 0;
1560  for (j = 0; j < n - 1; j++) {
1561  buff.fPols[indx++] = c + i;
1562  buff.fPols[indx++] = 4;
1563  buff.fPols[indx++] = nz * 2 * m + (2 * k + i * 1 + 2) * n + j + 1;
1564  buff.fPols[indx++] = (2 * k + i * 1 + 2) * m + j;
1565  buff.fPols[indx++] = nz * 2 * m + (2 * k + i * 1 + 2) * n + j;
1566  buff.fPols[indx++] = (2 * k + i * 1) * m + j;
1567  }
1568  if (specialCase) {
1569  buff.fPols[indx++] = c + i;
1570  buff.fPols[indx++] = 4;
1571  buff.fPols[indx++] = nz * 2 * m + (2 * k + i * 1 + 2) * n;
1572  buff.fPols[indx++] = (2 * k + i * 1 + 2) * m + j;
1573  buff.fPols[indx++] = nz * 2 * m + (2 * k + i * 1 + 2) * n + j;
1574  buff.fPols[indx++] = (2 * k + i * 1) * m + j;
1575  }
1576  i = 1;
1577  for (j = 0; j < n - 1; j++) {
1578  buff.fPols[indx++] = c + i;
1579  buff.fPols[indx++] = 4;
1580  buff.fPols[indx++] = (2 * k + i * 1) * m + j;
1581  buff.fPols[indx++] = nz * 2 * m + (2 * k + i * 1 + 2) * n + j;
1582  buff.fPols[indx++] = (2 * k + i * 1 + 2) * m + j;
1583  buff.fPols[indx++] = nz * 2 * m + (2 * k + i * 1 + 2) * n + j + 1;
1584  }
1585  if (specialCase) {
1586  buff.fPols[indx++] = c + i;
1587  buff.fPols[indx++] = 4;
1588  buff.fPols[indx++] = (2 * k + i * 1) * m + j;
1589  buff.fPols[indx++] = nz * 2 * m + (2 * k + i * 1 + 2) * n + j;
1590  buff.fPols[indx++] = (2 * k + i * 1 + 2) * m + j;
1591  buff.fPols[indx++] = nz * 2 * m + (2 * k + i * 1 + 2) * n;
1592  }
1593  }
1594 
1595  // left & right sections, number of polygons: 2*(nz-1)
1596  // special case number of polygons: 0
1597  if (!specialCase) {
1598  indx2 = nz * 2 * (n - 1);
1599  for (k = 0; k < (nz - 1); k++) {
1600  buff.fPols[indx++] = c + 2;
1601  buff.fPols[indx++] = 4;
1602  buff.fPols[indx++] = k == 0 ? indx2 : indx2 + 2 * nz * n + 2 * (k - 1);
1603  buff.fPols[indx++] = indx2 + 2 * (k + 1) * n;
1604  buff.fPols[indx++] = indx2 + 2 * nz * n + 2 * k;
1605  buff.fPols[indx++] = indx2 + (2 * k + 3) * n;
1606 
1607  buff.fPols[indx++] = c + 2;
1608  buff.fPols[indx++] = 4;
1609  buff.fPols[indx++] = k == 0 ? indx2 + n - 1 : indx2 + 2 * nz * n + 2 * (k - 1) + 1; // a
1610  buff.fPols[indx++] = indx2 + (2 * k + 3) * n + n - 1; // d
1611  buff.fPols[indx++] = indx2 + 2 * nz * n + 2 * k + 1; // c
1612  buff.fPols[indx++] = indx2 + 2 * (k + 1) * n + n - 1; // b
1613  }
1614  buff.fPols[indx - 8] = indx2 + n;
1615  buff.fPols[indx - 2] = indx2 + 2 * n - 1;
1616  }
1617 }
1618 
1619 ////////////////////////////////////////////////////////////////////////////////
1620 /// Computes projected pgon radius (inner or outer) corresponding to a given Z
1621 /// value. Fills corresponding coefficients of:
1622 /// `Rpg(z) = a + b*z`
1623 ///
1624 /// Note: ipl must be in range [0,fNz-2]
1625 
1626 Double_t TGeoPgon::Rpg(Double_t z, Int_t ipl, Bool_t inner, Double_t &a, Double_t &b) const
1627 {
1628  Double_t rpg;
1629  if (ipl < 0 || ipl > fNz - 2) {
1630  Fatal("Rpg", "Plane index parameter ipl=%i out of range\n", ipl);
1631  return 0;
1632  }
1633  Double_t dz = fZ[ipl + 1] - fZ[ipl];
1634  if (dz < TGeoShape::Tolerance()) {
1635  // radius-changing region
1636  rpg = (inner) ? TMath::Min(fRmin[ipl], fRmin[ipl + 1]) : TMath::Max(fRmax[ipl], fRmax[ipl + 1]);
1637  a = rpg;
1638  b = 0.;
1639  return rpg;
1640  }
1641  Double_t r1 = 0, r2 = 0;
1642  if (inner) {
1643  r1 = fRmin[ipl];
1644  r2 = fRmin[ipl + 1];
1645  } else {
1646  r1 = fRmax[ipl];
1647  r2 = fRmax[ipl + 1];
1648  }
1649  Double_t dzinv = 1. / dz;
1650  a = (r1 * fZ[ipl + 1] - r2 * fZ[ipl]) * dzinv;
1651  b = (r2 - r1) * dzinv;
1652  return (a + b * z);
1653 }
1654 
1655 ////////////////////////////////////////////////////////////////////////////////
1656 /// Computes projected distance at a given Z for a given ray inside a given sector
1657 /// and fills coefficients:
1658 /// `Rproj = a + b*z`
1659 
1660 Double_t TGeoPgon::Rproj(Double_t z, const Double_t *point, const Double_t *dir, Double_t cphi, Double_t sphi,
1661  Double_t &a, Double_t &b) const
1662 {
1663  if (TMath::Abs(dir[2]) < TGeoShape::Tolerance()) {
1664  a = b = TGeoShape::Big();
1665  return TGeoShape::Big();
1666  }
1667  Double_t invdirz = 1. / dir[2];
1668  a = ((point[0] * dir[2] - point[2] * dir[0]) * cphi + (point[1] * dir[2] - point[2] * dir[1]) * sphi) * invdirz;
1669  b = (dir[0] * cphi + dir[1] * sphi) * invdirz;
1670  return (a + b * z);
1671 }
1672 
1673 ////////////////////////////////////////////////////////////////////////////////
1674 /// Compute safety from POINT to segment between planes ipl, ipl+1 within safmin.
1675 
1676 Double_t TGeoPgon::SafetyToSegment(const Double_t *point, Int_t ipl, Int_t iphi, Bool_t in, Double_t safphi,
1677  Double_t safmin) const
1678 {
1679  Double_t saf[3];
1680  Double_t safe;
1681  Int_t i;
1682  Double_t r, rpgon, ta, calf;
1683  if (ipl < 0 || ipl > fNz - 2) return (safmin + 1.); // error in input plane
1684  // Get info about segment.
1685  Double_t dz = fZ[ipl + 1] - fZ[ipl];
1686  if (dz < 1E-9) return 1E9; // skip radius-changing segment
1687  Double_t znew = point[2] - 0.5 * (fZ[ipl] + fZ[ipl + 1]);
1688  saf[0] = 0.5 * dz - TMath::Abs(znew);
1689  if (-saf[0] > safmin) return TGeoShape::Big(); // means: stop checking further segments
1690  Double_t rmin1 = fRmin[ipl];
1691  Double_t rmax1 = fRmax[ipl];
1692  Double_t rmin2 = fRmin[ipl + 1];
1693  Double_t rmax2 = fRmax[ipl + 1];
1694  Double_t divphi = fDphi / fNedges;
1695  if (iphi < 0) {
1696  Double_t f = 1. / TMath::Cos(0.5 * divphi * TMath::DegToRad());
1697  rmax1 *= f;
1698  rmax2 *= f;
1699  r = TMath::Sqrt(point[0] * point[0] + point[1] * point[1]);
1700  Double_t ro1 = 0.5 * (rmin1 + rmin2);
1701  Double_t tg1 = (rmin2 - rmin1) / dz;
1702  Double_t cr1 = 1. / TMath::Sqrt(1. + tg1 * tg1);
1703  Double_t ro2 = 0.5 * (rmax1 + rmax2);
1704  Double_t tg2 = (rmax2 - rmax1) / dz;
1705  Double_t cr2 = 1. / TMath::Sqrt(1. + tg2 * tg2);
1706  Double_t rin = tg1 * znew + ro1;
1707  Double_t rout = tg2 * znew + ro2;
1708  saf[1] = (ro1 > 0) ? ((r - rin) * cr1) : TGeoShape::Big();
1709  saf[2] = (rout - r) * cr2;
1710  for (i = 0; i < 3; i++) saf[i] = -saf[i];
1711  safe = saf[TMath::LocMax(3, saf)];
1712  safe = TMath::Max(safe, safphi);
1713  if (safe < 0) safe = 0;
1714  return safe;
1715  }
1716  Double_t ph0 = (fPhi1 + divphi * (iphi + 0.5)) * TMath::DegToRad();
1717  r = point[0] * TMath::Cos(ph0) + point[1] * TMath::Sin(ph0);
1718  if (rmin1 + rmin2 > 1E-10) {
1719  ta = (rmin2 - rmin1) / dz;
1720  calf = 1. / TMath::Sqrt(1 + ta * ta);
1721  rpgon = rmin1 + (point[2] - fZ[ipl]) * ta;
1722  saf[1] = (r - rpgon) * calf;
1723  } else {
1724  saf[1] = TGeoShape::Big();
1725  }
1726  ta = (rmax2 - rmax1) / dz;
1727  calf = 1. / TMath::Sqrt(1 + ta * ta);
1728  rpgon = rmax1 + (point[2] - fZ[ipl]) * ta;
1729  saf[2] = (rpgon - r) * calf;
1730  if (in) {
1731  safe = saf[TMath::LocMin(3, saf)];
1732  safe = TMath::Min(safe, safphi);
1733  } else {
1734  for (i = 0; i < 3; i++) saf[i] = -saf[i];
1735  safe = saf[TMath::LocMax(3, saf)];
1736  safe = TMath::Max(safe, safphi);
1737  }
1738  if (safe < 0) safe = 0;
1739  return safe;
1740 }
1741 
1742 ////////////////////////////////////////////////////////////////////////////////
1743 /// computes the closest distance from given point to this shape, according
1744 /// to option. The matching point on the shape is stored in spoint.
1745 
1746 Double_t TGeoPgon::Safety(const Double_t *point, Bool_t in) const
1747 {
1748  Double_t safmin, saftmp, safphi;
1749  Double_t dz;
1750  Int_t ipl, iplane, iphi;
1751  LocatePhi(point, iphi);
1752  safphi = TGeoShape::SafetyPhi(point, in, fPhi1, fPhi1 + fDphi);
1753  if (in) {
1754  //---> point is inside pgon
1755  ipl = TMath::BinarySearch(fNz, fZ, point[2]);
1756  if (ipl == (fNz - 1)) return 0; // point on last Z boundary
1757  if (ipl < 0) return 0; // point on first Z boundary
1758  dz = 0.5 * (fZ[ipl + 1] - fZ[ipl]);
1759  if (dz < 1E-8) return 0;
1760  // Check safety for current segment
1761  safmin = SafetyToSegment(point, ipl, iphi, in, safphi);
1762  if (safmin > 1E10) {
1763  // something went wrong - point is not inside current segment
1764  return TGeoShape::Big();
1765  }
1766  if (safmin < 1E-6) return TMath::Abs(safmin); // point on radius-changing plane
1767  // check increasing iplanes
1768  iplane = ipl + 1;
1769  saftmp = 0.;
1770  while ((iplane < fNz - 1) && saftmp < 1E10) {
1771  saftmp = TMath::Abs(SafetyToSegment(point, iplane, iphi, kFALSE, safphi, safmin));
1772  if (saftmp < safmin) safmin = saftmp;
1773  iplane++;
1774  }
1775  // now decreasing nplanes
1776  iplane = ipl - 1;
1777  saftmp = 0.;
1778  while ((iplane >= 0) && saftmp < 1E10) {
1779  saftmp = TMath::Abs(SafetyToSegment(point, iplane, iphi, kFALSE, safphi, safmin));
1780  if (saftmp < safmin) safmin = saftmp;
1781  iplane--;
1782  }
1783  return safmin;
1784  }
1785  //---> point is outside pgon
1786  ipl = TMath::BinarySearch(fNz, fZ, point[2]);
1787  if (ipl < 0)
1788  ipl = 0;
1789  else if (ipl == fNz - 1)
1790  ipl = fNz - 2;
1791  dz = 0.5 * (fZ[ipl + 1] - fZ[ipl]);
1792  if (dz < 1E-8) {
1793  ipl++;
1794  if (ipl > fNz - 2) return 0.; // invalid last section
1795  dz = 0.5 * (fZ[ipl + 1] - fZ[ipl]);
1796  }
1797  // Check safety for current segment
1798  safmin = SafetyToSegment(point, ipl, iphi, kFALSE, safphi);
1799  if (safmin < 1E-6) return TMath::Abs(safmin); // point on radius-changing plane
1800  saftmp = 0.;
1801  // check increasing iplanes
1802  iplane = ipl + 1;
1803  saftmp = 0.;
1804  while ((iplane < fNz - 1) && saftmp < 1E10) {
1805  saftmp = TMath::Abs(SafetyToSegment(point, iplane, iphi, kFALSE, safphi, safmin));
1806  if (saftmp < safmin) safmin = saftmp;
1807  iplane++;
1808  }
1809  // now decreasing nplanes
1810  iplane = ipl - 1;
1811  saftmp = 0.;
1812  while ((iplane >= 0) && saftmp < 1E10) {
1813  saftmp = TMath::Abs(SafetyToSegment(point, iplane, iphi, kFALSE, safphi, safmin));
1814  if (saftmp < safmin) safmin = saftmp;
1815  iplane--;
1816  }
1817  return safmin;
1818 }
1819 
1820 ////////////////////////////////////////////////////////////////////////////////
1821 /// Save a primitive as a C++ statement(s) on output stream "out".
1822 
1823 void TGeoPgon::SavePrimitive(std::ostream &out, Option_t * /*option*/ /*= ""*/)
1824 {
1825  if (TObject::TestBit(kGeoSavePrimitive)) return;
1826  out << " // Shape: " << GetName() << " type: " << ClassName() << std::endl;
1827  out << " phi1 = " << fPhi1 << ";" << std::endl;
1828  out << " dphi = " << fDphi << ";" << std::endl;
1829  out << " nedges = " << fNedges << ";" << std::endl;
1830  out << " nz = " << fNz << ";" << std::endl;
1831  out << " TGeoPgon *pgon = new TGeoPgon(\"" << GetName() << "\",phi1,dphi,nedges,nz);" << std::endl;
1832  for (Int_t i = 0; i < fNz; i++) {
1833  out << " z = " << fZ[i] << ";" << std::endl;
1834  out << " rmin = " << fRmin[i] << ";" << std::endl;
1835  out << " rmax = " << fRmax[i] << ";" << std::endl;
1836  out << " pgon->DefineSection(" << i << ", z,rmin,rmax);" << std::endl;
1837  }
1838  out << " TGeoShape *" << GetPointerName() << " = pgon;" << std::endl;
1839  TObject::SetBit(TGeoShape::kGeoSavePrimitive);
1840 }
1841 
1842 ////////////////////////////////////////////////////////////////////////////////
1843 /// Set PGON dimensions starting from an array.
1844 
1845 void TGeoPgon::SetDimensions(Double_t *param)
1846 {
1847  fPhi1 = param[0];
1848  fDphi = param[1];
1849  fNedges = (Int_t)param[2];
1850  fNz = (Int_t)param[3];
1851  if (fNz < 2) {
1852  Error("SetDimensions", "Pgon %s: Number of Z sections must be > 2", GetName());
1853  return;
1854  }
1855  if (fRmin) delete[] fRmin;
1856  if (fRmax) delete[] fRmax;
1857  if (fZ) delete[] fZ;
1858  fRmin = new Double_t[fNz];
1859  fRmax = new Double_t[fNz];
1860  fZ = new Double_t[fNz];
1861  memset(fRmin, 0, fNz * sizeof(Double_t));
1862  memset(fRmax, 0, fNz * sizeof(Double_t));
1863  memset(fZ, 0, fNz * sizeof(Double_t));
1864  for (Int_t i = 0; i < fNz; i++) DefineSection(i, param[4 + 3 * i], param[5 + 3 * i], param[6 + 3 * i]);
1865 }
1866 
1867 ////////////////////////////////////////////////////////////////////////////////
1868 /// create polygone mesh points
1869 
1870 void TGeoPgon::SetPoints(Double_t *points) const
1871 {
1872  Double_t phi, dphi;
1873  Int_t n = fNedges + 1;
1874  dphi = fDphi / (n - 1);
1875  Double_t factor = 1. / TMath::Cos(TMath::DegToRad() * dphi / 2);
1876  Int_t i, j;
1877  Int_t indx = 0;
1878 
1879  if (points) {
1880  for (i = 0; i < fNz; i++) {
1881  for (j = 0; j < n; j++) {
1882  phi = (fPhi1 + j * dphi) * TMath::DegToRad();
1883  points[indx++] = factor * fRmin[i] * TMath::Cos(phi);
1884  points[indx++] = factor * fRmin[i] * TMath::Sin(phi);
1885  points[indx++] = fZ[i];
1886  }
1887  for (j = 0; j < n; j++) {
1888  phi = (fPhi1 + j * dphi) * TMath::DegToRad();
1889  points[indx++] = factor * fRmax[i] * TMath::Cos(phi);
1890  points[indx++] = factor * fRmax[i] * TMath::Sin(phi);
1891  points[indx++] = fZ[i];
1892  }
1893  }
1894  }
1895 }
1896 
1897 ////////////////////////////////////////////////////////////////////////////////
1898 /// create polygone mesh points
1899 
1900 void TGeoPgon::SetPoints(Float_t *points) const
1901 {
1902  Double_t phi, dphi;
1903  Int_t n = fNedges + 1;
1904  dphi = fDphi / (n - 1);
1905  Double_t factor = 1. / TMath::Cos(TMath::DegToRad() * dphi / 2);
1906  Int_t i, j;
1907  Int_t indx = 0;
1908 
1909  if (points) {
1910  for (i = 0; i < fNz; i++) {
1911  for (j = 0; j < n; j++) {
1912  phi = (fPhi1 + j * dphi) * TMath::DegToRad();
1913  points[indx++] = factor * fRmin[i] * TMath::Cos(phi);
1914  points[indx++] = factor * fRmin[i] * TMath::Sin(phi);
1915  points[indx++] = fZ[i];
1916  }
1917  for (j = 0; j < n; j++) {
1918  phi = (fPhi1 + j * dphi) * TMath::DegToRad();
1919  points[indx++] = factor * fRmax[i] * TMath::Cos(phi);
1920  points[indx++] = factor * fRmax[i] * TMath::Sin(phi);
1921  points[indx++] = fZ[i];
1922  }
1923  }
1924  }
1925 }
1926 
1927 ////////////////////////////////////////////////////////////////////////////////
1928 /// Returns numbers of vertices, segments and polygons composing the shape mesh.
1929 
1930 void TGeoPgon::GetMeshNumbers(Int_t &nvert, Int_t &nsegs, Int_t &npols) const
1931 {
1932  Int_t n = fNedges + 1;
1933  Int_t nz = GetNz();
1934  nvert = nz * 2 * n;
1935  Bool_t specialCase = (TGeoShape::IsSameWithinTolerance(GetDphi(), 360));
1936  nsegs = 4 * (nz * n - 1 + (specialCase == kTRUE));
1937  npols = 2 * (nz * n - 1 + (specialCase == kTRUE));
1938 }
1939 
1940 ////////////////////////////////////////////////////////////////////////////////
1941 /// Return number of vertices of the mesh representation
1942 
1943 Int_t TGeoPgon::GetNmeshVertices() const
1944 {
1945  Int_t n = fNedges + 1;
1946  Int_t numPoints = fNz * 2 * n;
1947  return numPoints;
1948 }
1949 
1950 ////////////////////////////////////////////////////////////////////////////////
1951 /// fill size of this 3-D object
1952 
1953 void TGeoPgon::Sizeof3D() const
1954 {
1955 }
1956 
1957 ////////////////////////////////////////////////////////////////////////////////
1958 /// Fills a static 3D buffer and returns a reference.
1959 
1960 const TBuffer3D &TGeoPgon::GetBuffer3D(Int_t reqSections, Bool_t localFrame) const
1961 {
1962  static TBuffer3D buffer(TBuffer3DTypes::kGeneric);
1963 
1964  TGeoBBox::FillBuffer3D(buffer, reqSections, localFrame);
1965 
1966  if (reqSections & TBuffer3D::kRawSizes) {
1967  const Int_t n = GetNsegments() + 1;
1968  Int_t nz = GetNz();
1969  Int_t nbPnts = nz * 2 * n;
1970  if (nz >= 2 && nbPnts > 0) {
1971  Bool_t specialCase = (TGeoShape::IsSameWithinTolerance(GetDphi(), 360));
1972  Int_t nbSegs = 4 * (nz * n - 1 + (specialCase == kTRUE));
1973  Int_t nbPols = 2 * (nz * n - 1 + (specialCase == kTRUE));
1974  if (buffer.SetRawSizes(nbPnts, 3 * nbPnts, nbSegs, 3 * nbSegs, nbPols, 6 * nbPols)) {
1975  buffer.SetSectionsValid(TBuffer3D::kRawSizes);
1976  }
1977  }
1978  }
1979  // TODO: Push down to TGeoShape?? Would have to do raw sizes set first..
1980  // can rest of TGeoShape be deferred until after this?
1981  if ((reqSections & TBuffer3D::kRaw) && buffer.SectionsValid(TBuffer3D::kRawSizes)) {
1982  SetPoints(buffer.fPnts);
1983  if (!buffer.fLocalFrame) {
1984  TransformPoints(buffer.fPnts, buffer.NbPnts());
1985  }
1986 
1987  SetSegsAndPols(buffer);
1988  buffer.SetSectionsValid(TBuffer3D::kRaw);
1989  }
1990 
1991  return buffer;
1992 }
1993 
1994 ////////////////////////////////////////////////////////////////////////////////
1995 /// Check the inside status for each of the points in the array.
1996 /// Input: Array of point coordinates + vector size
1997 /// Output: Array of Booleans for the inside of each point
1998 
1999 void TGeoPgon::Contains_v(const Double_t *points, Bool_t *inside, Int_t vecsize) const
2000 {
2001  for (Int_t i = 0; i < vecsize; i++) inside[i] = Contains(&points[3 * i]);
2002 }
2003 
2004 ////////////////////////////////////////////////////////////////////////////////
2005 /// Compute the normal for an array o points so that norm.dot.dir is positive
2006 /// Input: Arrays of point coordinates and directions + vector size
2007 /// Output: Array of normal directions
2008 
2009 void TGeoPgon::ComputeNormal_v(const Double_t *points, const Double_t *dirs, Double_t *norms, Int_t vecsize)
2010 {
2011  for (Int_t i = 0; i < vecsize; i++) ComputeNormal(&points[3 * i], &dirs[3 * i], &norms[3 * i]);
2012 }
2013 
2014 ////////////////////////////////////////////////////////////////////////////////
2015 /// Compute distance from array of input points having directions specified by dirs. Store output in dists
2016 
2017 void TGeoPgon::DistFromInside_v(const Double_t *points, const Double_t *dirs, Double_t *dists, Int_t vecsize,
2018  Double_t *step) const
2019 {
2020  for (Int_t i = 0; i < vecsize; i++) dists[i] = DistFromInside(&points[3 * i], &dirs[3 * i], 3, step[i]);
2021 }
2022 
2023 ////////////////////////////////////////////////////////////////////////////////
2024 /// Compute distance from array of input points having directions specified by dirs. Store output in dists
2025 
2026 void TGeoPgon::DistFromOutside_v(const Double_t *points, const Double_t *dirs, Double_t *dists, Int_t vecsize,
2027  Double_t *step) const
2028 {
2029  for (Int_t i = 0; i < vecsize; i++) dists[i] = DistFromOutside(&points[3 * i], &dirs[3 * i], 3, step[i]);
2030 }
2031 
2032 ////////////////////////////////////////////////////////////////////////////////
2033 /// Compute safe distance from each of the points in the input array.
2034 /// Input: Array of point coordinates, array of statuses for these points, size of the arrays
2035 /// Output: Safety values
2036 
2037 void TGeoPgon::Safety_v(const Double_t *points, const Bool_t *inside, Double_t *safe, Int_t vecsize) const
2038 {
2039  for (Int_t i = 0; i < vecsize; i++) safe[i] = Safety(&points[3 * i], inside[i]);
2040 }