Logo ROOT  
Reference Guide
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
16A 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
25Begin_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}
47End_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
64
65 ////////////////////////////////////////////////////////////////////////////////
66 /// Constructor.
67
69 : fIntBuffer(0), fDblBuffer(0)
70{
71}
72
73////////////////////////////////////////////////////////////////////////////////
74/// Destructor.
75
77{
78 delete[] fIntBuffer;
79 delete[] fDblBuffer;
80}
81
82////////////////////////////////////////////////////////////////////////////////
83
85{
87 return *fThreadData[tid];
88}
89
90////////////////////////////////////////////////////////////////////////////////
91
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
108{
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
126{
128 fNedges = 0;
129 fThreadSize = 0;
130}
131
132////////////////////////////////////////////////////////////////////////////////
133/// Default constructor
134
135TGeoPgon::TGeoPgon(Double_t phi, Double_t dphi, Int_t nedges, Int_t nz) : TGeoPcon(phi, dphi, nz)
136{
138 fNedges = nedges;
139 fThreadSize = 0;
141}
142
143////////////////////////////////////////////////////////////////////////////////
144/// Default constructor
145
146TGeoPgon::TGeoPgon(const char *name, Double_t phi, Double_t dphi, Int_t nedges, Int_t nz)
147 : TGeoPcon(name, phi, dphi, nz)
148{
150 fNedges = nedges;
151 fThreadSize = 0;
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
167{
169 SetDimensions(param);
170 ComputeBBox();
171 fThreadSize = 0;
173}
174
175////////////////////////////////////////////////////////////////////////////////
176/// destructor
177
179{
181}
182
183////////////////////////////////////////////////////////////////////////////////
184/// Computes capacity of the shape in [length^3]
185
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
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(),
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);
272}
273
274////////////////////////////////////////////////////////////////////////////////
275/// Compute normal to closest surface from POINT.
276
277void 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
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
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);
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);
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
519void 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
530Int_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
610Bool_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;
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
683Bool_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;
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
764Bool_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
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
900Bool_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
989Bool_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]);
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");
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;
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]);
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
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;
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
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);
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;
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
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
1299TGeoVolume *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
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];
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
1409{
1410 printf("*** Shape %s: TGeoPgon ***\n", GetName());
1411 printf(" Nedges = %i\n", fNedges);
1413}
1414
1415////////////////////////////////////////////////////////////////////////////////
1416/// Creates a TBuffer3D describing *this* shape.
1417/// Coordinates are in local reference frame.
1418
1420{
1421 Int_t nbPnts, nbSegs, nbPols;
1422 GetMeshNumbers(nbPnts, nbSegs, nbPols);
1423
1424 if (nbPnts <= 0)
1425 return nullptr;
1426
1427 TBuffer3D *buff =
1428 new TBuffer3D(TBuffer3DTypes::kGeneric, nbPnts, 3 * nbPnts, nbSegs, 3 * nbSegs, nbPols, 6 * nbPols);
1429 if (buff) {
1430 SetPoints(buff->fPnts);
1431 SetSegsAndPols(*buff);
1432 }
1433
1434 return buff;
1435}
1436
1437////////////////////////////////////////////////////////////////////////////////
1438/// Fill TBuffer3D structure for segments and polygons.
1439
1441{
1442 if (!HasInsideSurface()) {
1444 return;
1445 }
1446
1447 Int_t i, j;
1448 const Int_t n = GetNedges() + 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 = TGeoShape::IsSameWithinTolerance(dphi, 360);
1455
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 ? 1 : 0);
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////////////////////////////////////////////////////////////////////////////////
1621/// Fill TBuffer3D structure for segments and polygons, when no inner surface exists
1622
1624{
1625 const Int_t n = GetNedges() + 1;
1626 const Int_t nz = GetNz();
1627 const Int_t nbPnts = nz * n + 2;
1628
1629 if ((nz < 2) || (nbPnts <= 0) || (n < 2)) return;
1630
1631 Int_t c = GetBasicColor();
1632
1633 Int_t indx = 0, indx1 = 0, indx2 = 0, i, j;
1634
1635 // outside circles, number of segments: nz*n
1636 for (i = 0; i < nz; i++) {
1637 indx2 = i * n;
1638 for (j = 1; j < n; j++) {
1639 buff.fSegs[indx++] = c;
1640 buff.fSegs[indx++] = indx2 + j - 1;
1641 buff.fSegs[indx++] = indx2 + j % (n-1);
1642 }
1643 }
1644
1645 indx2 = 0;
1646 // bottom lines
1647 for (j = 0; j < n; j++) {
1648 buff.fSegs[indx++] = c;
1649 buff.fSegs[indx++] = indx2 + j % (n-1);
1650 buff.fSegs[indx++] = nbPnts - 2;
1651 }
1652
1653 indx2 = (nz-1) * n;
1654 // top lines
1655 for (j = 0; j < n; j++) {
1656 buff.fSegs[indx++] = c;
1657 buff.fSegs[indx++] = indx2 + j % (n-1);
1658 buff.fSegs[indx++] = nbPnts - 1;
1659 }
1660
1661 // outside cylinders, number of segments: (nz-1)*n
1662 for (i = 0; i < (nz - 1); i++) {
1663 // outside cylinder
1664 indx2 = i * n;
1665 for (j = 0; j < n; j++) {
1666 buff.fSegs[indx++] = c;
1667 buff.fSegs[indx++] = indx2 + j % (n-1);
1668 buff.fSegs[indx++] = indx2 + n + j % (n-1);
1669 }
1670 }
1671
1672 indx = 0;
1673
1674 // bottom cap
1675 indx1 = 0; // start of first z layer
1676 indx2 = nz*(n-1);
1677 for (j = 0; j < n - 1; j++) {
1678 buff.fPols[indx++] = c;
1679 buff.fPols[indx++] = 3;
1680 buff.fPols[indx++] = indx1 + j;
1681 buff.fPols[indx++] = indx2 + (j+1)%(n-1);
1682 buff.fPols[indx++] = indx2 + j;
1683 }
1684
1685 // top cap
1686 indx1 = (nz-1)*(n-1); // start last z layer
1687 indx2 = nz*(n-1) + n;
1688 for (j = 0; j < n - 1; j++) {
1689 buff.fPols[indx++] = c;
1690 buff.fPols[indx++] = 3;
1691 buff.fPols[indx++] = indx1 + j; // last z layer
1692 buff.fPols[indx++] = indx2 + j;
1693 buff.fPols[indx++] = indx2 + (j+1)%(n-1);
1694 }
1695
1696 // outside, number of polygons: (nz-1)*(n-1)
1697 for (Int_t k = 0; k < (nz - 1); k++) {
1698 indx1 = k*(n-1);
1699 indx2 = nz*(n-1) + n*2 + k*n;
1700 for (j = 0; j < n-1; j++) {
1701 buff.fPols[indx++] = c;
1702 buff.fPols[indx++] = 4;
1703 buff.fPols[indx++] = indx1 + j;
1704 buff.fPols[indx++] = indx2 + j;
1705 buff.fPols[indx++] = indx1 + j + (n-1);
1706 buff.fPols[indx++] = indx2 + (j+1)%(n-1);
1707 }
1708 }
1709}
1710
1711
1712////////////////////////////////////////////////////////////////////////////////
1713/// Computes projected pgon radius (inner or outer) corresponding to a given Z
1714/// value. Fills corresponding coefficients of:
1715/// `Rpg(z) = a + b*z`
1716///
1717/// Note: ipl must be in range [0,fNz-2]
1718
1720{
1721 Double_t rpg;
1722 if (ipl < 0 || ipl > fNz - 2) {
1723 Fatal("Rpg", "Plane index parameter ipl=%i out of range\n", ipl);
1724 return 0;
1725 }
1726 Double_t dz = fZ[ipl + 1] - fZ[ipl];
1727 if (dz < TGeoShape::Tolerance()) {
1728 // radius-changing region
1729 rpg = (inner) ? TMath::Min(fRmin[ipl], fRmin[ipl + 1]) : TMath::Max(fRmax[ipl], fRmax[ipl + 1]);
1730 a = rpg;
1731 b = 0.;
1732 return rpg;
1733 }
1734 Double_t r1 = 0, r2 = 0;
1735 if (inner) {
1736 r1 = fRmin[ipl];
1737 r2 = fRmin[ipl + 1];
1738 } else {
1739 r1 = fRmax[ipl];
1740 r2 = fRmax[ipl + 1];
1741 }
1742 Double_t dzinv = 1. / dz;
1743 a = (r1 * fZ[ipl + 1] - r2 * fZ[ipl]) * dzinv;
1744 b = (r2 - r1) * dzinv;
1745 return (a + b * z);
1746}
1747
1748////////////////////////////////////////////////////////////////////////////////
1749/// Computes projected distance at a given Z for a given ray inside a given sector
1750/// and fills coefficients:
1751/// `Rproj = a + b*z`
1752
1753Double_t TGeoPgon::Rproj(Double_t z, const Double_t *point, const Double_t *dir, Double_t cphi, Double_t sphi,
1754 Double_t &a, Double_t &b) const
1755{
1756 if (TMath::Abs(dir[2]) < TGeoShape::Tolerance()) {
1757 a = b = TGeoShape::Big();
1758 return TGeoShape::Big();
1759 }
1760 Double_t invdirz = 1. / dir[2];
1761 a = ((point[0] * dir[2] - point[2] * dir[0]) * cphi + (point[1] * dir[2] - point[2] * dir[1]) * sphi) * invdirz;
1762 b = (dir[0] * cphi + dir[1] * sphi) * invdirz;
1763 return (a + b * z);
1764}
1765
1766////////////////////////////////////////////////////////////////////////////////
1767/// Compute safety from POINT to segment between planes ipl, ipl+1 within safmin.
1768
1770 Double_t safmin) const
1771{
1772 Double_t saf[3];
1773 Double_t safe;
1774 Int_t i;
1775 Double_t r, rpgon, ta, calf;
1776 if (ipl < 0 || ipl > fNz - 2) return (safmin + 1.); // error in input plane
1777 // Get info about segment.
1778 Double_t dz = fZ[ipl + 1] - fZ[ipl];
1779 if (dz < 1E-9) return 1E9; // skip radius-changing segment
1780 Double_t znew = point[2] - 0.5 * (fZ[ipl] + fZ[ipl + 1]);
1781 saf[0] = 0.5 * dz - TMath::Abs(znew);
1782 if (-saf[0] > safmin) return TGeoShape::Big(); // means: stop checking further segments
1783 Double_t rmin1 = fRmin[ipl];
1784 Double_t rmax1 = fRmax[ipl];
1785 Double_t rmin2 = fRmin[ipl + 1];
1786 Double_t rmax2 = fRmax[ipl + 1];
1787 Double_t divphi = fDphi / fNedges;
1788 if (iphi < 0) {
1789 Double_t f = 1. / TMath::Cos(0.5 * divphi * TMath::DegToRad());
1790 rmax1 *= f;
1791 rmax2 *= f;
1792 r = TMath::Sqrt(point[0] * point[0] + point[1] * point[1]);
1793 Double_t ro1 = 0.5 * (rmin1 + rmin2);
1794 Double_t tg1 = (rmin2 - rmin1) / dz;
1795 Double_t cr1 = 1. / TMath::Sqrt(1. + tg1 * tg1);
1796 Double_t ro2 = 0.5 * (rmax1 + rmax2);
1797 Double_t tg2 = (rmax2 - rmax1) / dz;
1798 Double_t cr2 = 1. / TMath::Sqrt(1. + tg2 * tg2);
1799 Double_t rin = tg1 * znew + ro1;
1800 Double_t rout = tg2 * znew + ro2;
1801 saf[1] = (ro1 > 0) ? ((r - rin) * cr1) : TGeoShape::Big();
1802 saf[2] = (rout - r) * cr2;
1803 for (i = 0; i < 3; i++) saf[i] = -saf[i];
1804 safe = saf[TMath::LocMax(3, saf)];
1805 safe = TMath::Max(safe, safphi);
1806 if (safe < 0) safe = 0;
1807 return safe;
1808 }
1809 Double_t ph0 = (fPhi1 + divphi * (iphi + 0.5)) * TMath::DegToRad();
1810 r = point[0] * TMath::Cos(ph0) + point[1] * TMath::Sin(ph0);
1811 if (rmin1 + rmin2 > 1E-10) {
1812 ta = (rmin2 - rmin1) / dz;
1813 calf = 1. / TMath::Sqrt(1 + ta * ta);
1814 rpgon = rmin1 + (point[2] - fZ[ipl]) * ta;
1815 saf[1] = (r - rpgon) * calf;
1816 } else {
1817 saf[1] = TGeoShape::Big();
1818 }
1819 ta = (rmax2 - rmax1) / dz;
1820 calf = 1. / TMath::Sqrt(1 + ta * ta);
1821 rpgon = rmax1 + (point[2] - fZ[ipl]) * ta;
1822 saf[2] = (rpgon - r) * calf;
1823 if (in) {
1824 safe = saf[TMath::LocMin(3, saf)];
1825 safe = TMath::Min(safe, safphi);
1826 } else {
1827 for (i = 0; i < 3; i++) saf[i] = -saf[i];
1828 safe = saf[TMath::LocMax(3, saf)];
1829 safe = TMath::Max(safe, safphi);
1830 }
1831 if (safe < 0) safe = 0;
1832 return safe;
1833}
1834
1835////////////////////////////////////////////////////////////////////////////////
1836/// computes the closest distance from given point to this shape, according
1837/// to option. The matching point on the shape is stored in spoint.
1838
1840{
1841 Double_t safmin, saftmp, safphi;
1842 Double_t dz;
1843 Int_t ipl, iplane, iphi;
1844 LocatePhi(point, iphi);
1845 safphi = TGeoShape::SafetyPhi(point, in, fPhi1, fPhi1 + fDphi);
1846 if (in) {
1847 //---> point is inside pgon
1848 ipl = TMath::BinarySearch(fNz, fZ, point[2]);
1849 if (ipl == (fNz - 1)) return 0; // point on last Z boundary
1850 if (ipl < 0) return 0; // point on first Z boundary
1851 dz = 0.5 * (fZ[ipl + 1] - fZ[ipl]);
1852 if (dz < 1E-8) return 0;
1853 // Check safety for current segment
1854 safmin = SafetyToSegment(point, ipl, iphi, in, safphi);
1855 if (safmin > 1E10) {
1856 // something went wrong - point is not inside current segment
1857 return TGeoShape::Big();
1858 }
1859 if (safmin < 1E-6) return TMath::Abs(safmin); // point on radius-changing plane
1860 // check increasing iplanes
1861 iplane = ipl + 1;
1862 saftmp = 0.;
1863 while ((iplane < fNz - 1) && saftmp < 1E10) {
1864 saftmp = TMath::Abs(SafetyToSegment(point, iplane, iphi, kFALSE, safphi, safmin));
1865 if (saftmp < safmin) safmin = saftmp;
1866 iplane++;
1867 }
1868 // now decreasing nplanes
1869 iplane = ipl - 1;
1870 saftmp = 0.;
1871 while ((iplane >= 0) && saftmp < 1E10) {
1872 saftmp = TMath::Abs(SafetyToSegment(point, iplane, iphi, kFALSE, safphi, safmin));
1873 if (saftmp < safmin) safmin = saftmp;
1874 iplane--;
1875 }
1876 return safmin;
1877 }
1878 //---> point is outside pgon
1879 ipl = TMath::BinarySearch(fNz, fZ, point[2]);
1880 if (ipl < 0)
1881 ipl = 0;
1882 else if (ipl == fNz - 1)
1883 ipl = fNz - 2;
1884 dz = 0.5 * (fZ[ipl + 1] - fZ[ipl]);
1885 if (dz < 1E-8) {
1886 ipl++;
1887 if (ipl > fNz - 2) return 0.; // invalid last section
1888 dz = 0.5 * (fZ[ipl + 1] - fZ[ipl]);
1889 }
1890 // Check safety for current segment
1891 safmin = SafetyToSegment(point, ipl, iphi, kFALSE, safphi);
1892 if (safmin < 1E-6) return TMath::Abs(safmin); // point on radius-changing plane
1893 saftmp = 0.;
1894 // check increasing iplanes
1895 iplane = ipl + 1;
1896 saftmp = 0.;
1897 while ((iplane < fNz - 1) && saftmp < 1E10) {
1898 saftmp = TMath::Abs(SafetyToSegment(point, iplane, iphi, kFALSE, safphi, safmin));
1899 if (saftmp < safmin) safmin = saftmp;
1900 iplane++;
1901 }
1902 // now decreasing nplanes
1903 iplane = ipl - 1;
1904 saftmp = 0.;
1905 while ((iplane >= 0) && saftmp < 1E10) {
1906 saftmp = TMath::Abs(SafetyToSegment(point, iplane, iphi, kFALSE, safphi, safmin));
1907 if (saftmp < safmin) safmin = saftmp;
1908 iplane--;
1909 }
1910 return safmin;
1911}
1912
1913////////////////////////////////////////////////////////////////////////////////
1914/// Save a primitive as a C++ statement(s) on output stream "out".
1915
1916void TGeoPgon::SavePrimitive(std::ostream &out, Option_t * /*option*/ /*= ""*/)
1917{
1919 out << " // Shape: " << GetName() << " type: " << ClassName() << std::endl;
1920 out << " phi1 = " << fPhi1 << ";" << std::endl;
1921 out << " dphi = " << fDphi << ";" << std::endl;
1922 out << " nedges = " << fNedges << ";" << std::endl;
1923 out << " nz = " << fNz << ";" << std::endl;
1924 out << " TGeoPgon *pgon = new TGeoPgon(\"" << GetName() << "\",phi1,dphi,nedges,nz);" << std::endl;
1925 for (Int_t i = 0; i < fNz; i++) {
1926 out << " z = " << fZ[i] << ";" << std::endl;
1927 out << " rmin = " << fRmin[i] << ";" << std::endl;
1928 out << " rmax = " << fRmax[i] << ";" << std::endl;
1929 out << " pgon->DefineSection(" << i << ", z,rmin,rmax);" << std::endl;
1930 }
1931 out << " TGeoShape *" << GetPointerName() << " = pgon;" << std::endl;
1933}
1934
1935////////////////////////////////////////////////////////////////////////////////
1936/// Set PGON dimensions starting from an array.
1937
1939{
1940 fPhi1 = param[0];
1941 fDphi = param[1];
1942 fNedges = (Int_t)param[2];
1943 fNz = (Int_t)param[3];
1944 if (fNz < 2) {
1945 Error("SetDimensions", "Pgon %s: Number of Z sections must be > 2", GetName());
1946 return;
1947 }
1948 if (fRmin) delete[] fRmin;
1949 if (fRmax) delete[] fRmax;
1950 if (fZ) delete[] fZ;
1951 fRmin = new Double_t[fNz];
1952 fRmax = new Double_t[fNz];
1953 fZ = new Double_t[fNz];
1954 memset(fRmin, 0, fNz * sizeof(Double_t));
1955 memset(fRmax, 0, fNz * sizeof(Double_t));
1956 memset(fZ, 0, fNz * sizeof(Double_t));
1957 for (Int_t i = 0; i < fNz; i++) DefineSection(i, param[4 + 3 * i], param[5 + 3 * i], param[6 + 3 * i]);
1958}
1959
1960////////////////////////////////////////////////////////////////////////////////
1961/// create polygone mesh points
1962
1964{
1965 Double_t phi, dphi;
1966 Int_t n = fNedges + 1;
1967 dphi = fDphi / (n - 1);
1968 Double_t factor = 1. / TMath::Cos(TMath::DegToRad() * dphi / 2);
1969 Int_t i, j;
1970 Int_t indx = 0;
1971
1972 Bool_t hasInside = HasInsideSurface();
1973
1974 if (points) {
1975 for (i = 0; i < GetNz(); i++) {
1976 if (hasInside)
1977 for (j = 0; j < n; j++) {
1978 phi = (fPhi1 + j * dphi) * TMath::DegToRad();
1979 points[indx++] = factor * fRmin[i] * TMath::Cos(phi);
1980 points[indx++] = factor * fRmin[i] * TMath::Sin(phi);
1981 points[indx++] = fZ[i];
1982 }
1983 for (j = 0; j < n; j++) {
1984 phi = (fPhi1 + j * dphi) * TMath::DegToRad();
1985 points[indx++] = factor * fRmax[i] * TMath::Cos(phi);
1986 points[indx++] = factor * fRmax[i] * TMath::Sin(phi);
1987 points[indx++] = fZ[i];
1988 }
1989 }
1990
1991 if (!hasInside) {
1992 points[indx++] = 0;
1993 points[indx++] = 0;
1994 points[indx++] = fZ[0];
1995
1996 points[indx++] = 0;
1997 points[indx++] = 0;
1998 points[indx++] = fZ[GetNz()-1];
1999 }
2000 }
2001}
2002
2003////////////////////////////////////////////////////////////////////////////////
2004/// create polygone mesh points
2005
2007{
2008 Double_t phi, dphi;
2009 Int_t n = fNedges + 1;
2010 dphi = fDphi / (n - 1);
2011 Double_t factor = 1. / TMath::Cos(TMath::DegToRad() * dphi / 2);
2012 Int_t i, j;
2013 Int_t indx = 0;
2014
2015 Bool_t hasInside = HasInsideSurface();
2016
2017 if (points) {
2018 for (i = 0; i < fNz; i++) {
2019 if (hasInside)
2020 for (j = 0; j < n; j++) {
2021 phi = (fPhi1 + j * dphi) * TMath::DegToRad();
2022 points[indx++] = factor * fRmin[i] * TMath::Cos(phi);
2023 points[indx++] = factor * fRmin[i] * TMath::Sin(phi);
2024 points[indx++] = fZ[i];
2025 }
2026 for (j = 0; j < n; j++) {
2027 phi = (fPhi1 + j * dphi) * TMath::DegToRad();
2028 points[indx++] = factor * fRmax[i] * TMath::Cos(phi);
2029 points[indx++] = factor * fRmax[i] * TMath::Sin(phi);
2030 points[indx++] = fZ[i];
2031 }
2032 }
2033
2034 if (!hasInside) {
2035 points[indx++] = 0;
2036 points[indx++] = 0;
2037 points[indx++] = fZ[0];
2038
2039 points[indx++] = 0;
2040 points[indx++] = 0;
2041 points[indx++] = fZ[GetNz()-1];
2042 }
2043 }
2044}
2045
2046////////////////////////////////////////////////////////////////////////////////
2047/// Returns numbers of vertices, segments and polygons composing the shape mesh.
2048
2049void TGeoPgon::GetMeshNumbers(Int_t &nvert, Int_t &nsegs, Int_t &npols) const
2050{
2051 nvert = nsegs = npols = 0;
2052
2053 Int_t n = GetNedges() + 1;
2054 Int_t nz = GetNz();
2055
2056 if (nz < 2) return;
2057
2058 if (HasInsideSurface()) {
2059 Bool_t specialCase = TGeoShape::IsSameWithinTolerance(GetDphi(), 360);
2060 nvert = nz * 2 * n;
2061 nsegs = 4 * (nz * n - 1 + (specialCase ? 1 : 0));
2062 npols = 2 * (nz * n - 1 + (specialCase ? 1 : 0));
2063 } else {
2064 nvert = nz * n + 2;
2065 nsegs = nz * (n - 1) + n * 2 + (nz - 1) * n;
2066 npols = 2 * (n - 1) + (nz - 1) * (n - 1);
2067 }
2068}
2069
2070////////////////////////////////////////////////////////////////////////////////
2071/// Return number of vertices of the mesh representation
2072
2074{
2075 Int_t nvert, nsegs, npols;
2076
2077 GetMeshNumbers(nvert, nsegs, npols);
2078
2079 return nvert;
2080}
2081
2082////////////////////////////////////////////////////////////////////////////////
2083/// fill size of this 3-D object
2084
2086{
2087}
2088
2089////////////////////////////////////////////////////////////////////////////////
2090/// Fills a static 3D buffer and returns a reference.
2091
2092const TBuffer3D &TGeoPgon::GetBuffer3D(Int_t reqSections, Bool_t localFrame) const
2093{
2094 static TBuffer3D buffer(TBuffer3DTypes::kGeneric);
2095
2096 TGeoBBox::FillBuffer3D(buffer, reqSections, localFrame);
2097
2098 if (reqSections & TBuffer3D::kRawSizes) {
2099 Int_t nbPnts, nbSegs, nbPols;
2100 GetMeshNumbers(nbPnts, nbSegs, nbPols);
2101 if (nbPnts > 0) {
2102 if (buffer.SetRawSizes(nbPnts, 3 * nbPnts, nbSegs, 3 * nbSegs, nbPols, 6 * nbPols)) {
2104 }
2105 }
2106 }
2107 // TODO: Push down to TGeoShape?? Would have to do raw sizes set first..
2108 // can rest of TGeoShape be deferred until after this?
2109 if ((reqSections & TBuffer3D::kRaw) && buffer.SectionsValid(TBuffer3D::kRawSizes)) {
2110 SetPoints(buffer.fPnts);
2111 if (!buffer.fLocalFrame) {
2112 TransformPoints(buffer.fPnts, buffer.NbPnts());
2113 }
2114
2115 SetSegsAndPols(buffer);
2117 }
2118
2119 return buffer;
2120}
2121
2122////////////////////////////////////////////////////////////////////////////////
2123/// Check the inside status for each of the points in the array.
2124/// Input: Array of point coordinates + vector size
2125/// Output: Array of Booleans for the inside of each point
2126
2127void TGeoPgon::Contains_v(const Double_t *points, Bool_t *inside, Int_t vecsize) const
2128{
2129 for (Int_t i = 0; i < vecsize; i++) inside[i] = Contains(&points[3 * i]);
2130}
2131
2132////////////////////////////////////////////////////////////////////////////////
2133/// Compute the normal for an array o points so that norm.dot.dir is positive
2134/// Input: Arrays of point coordinates and directions + vector size
2135/// Output: Array of normal directions
2136
2137void TGeoPgon::ComputeNormal_v(const Double_t *points, const Double_t *dirs, Double_t *norms, Int_t vecsize)
2138{
2139 for (Int_t i = 0; i < vecsize; i++) ComputeNormal(&points[3 * i], &dirs[3 * i], &norms[3 * i]);
2140}
2141
2142////////////////////////////////////////////////////////////////////////////////
2143/// Compute distance from array of input points having directions specified by dirs. Store output in dists
2144
2145void TGeoPgon::DistFromInside_v(const Double_t *points, const Double_t *dirs, Double_t *dists, Int_t vecsize,
2146 Double_t *step) const
2147{
2148 for (Int_t i = 0; i < vecsize; i++) dists[i] = DistFromInside(&points[3 * i], &dirs[3 * i], 3, step[i]);
2149}
2150
2151////////////////////////////////////////////////////////////////////////////////
2152/// Compute distance from array of input points having directions specified by dirs. Store output in dists
2153
2154void TGeoPgon::DistFromOutside_v(const Double_t *points, const Double_t *dirs, Double_t *dists, Int_t vecsize,
2155 Double_t *step) const
2156{
2157 for (Int_t i = 0; i < vecsize; i++) dists[i] = DistFromOutside(&points[3 * i], &dirs[3 * i], 3, step[i]);
2158}
2159
2160////////////////////////////////////////////////////////////////////////////////
2161/// Compute safe distance from each of the points in the input array.
2162/// Input: Array of point coordinates, array of statuses for these points, size of the arrays
2163/// Output: Safety values
2164
2165void TGeoPgon::Safety_v(const Double_t *points, const Bool_t *inside, Double_t *safe, Int_t vecsize) const
2166{
2167 for (Int_t i = 0; i < vecsize; i++) safe[i] = Safety(&points[3 * i], inside[i]);
2168}
ROOT::R::TRInterface & r
Definition: Object.C:4
#define b(i)
Definition: RSha256.hxx:100
#define f(i)
Definition: RSha256.hxx:104
#define c(i)
Definition: RSha256.hxx:101
#define s1(x)
Definition: RSha256.hxx:91
int Int_t
Definition: RtypesCore.h:43
const Bool_t kFALSE
Definition: RtypesCore.h:90
double Double_t
Definition: RtypesCore.h:57
float Float_t
Definition: RtypesCore.h:55
const Bool_t kTRUE
Definition: RtypesCore.h:89
const char Option_t
Definition: RtypesCore.h:64
#define ClassImp(name)
Definition: Rtypes.h:361
XFontStruct * id
Definition: TGX11.cxx:108
char name[80]
Definition: TGX11.cxx:109
R__EXTERN TGeoManager * gGeoManager
Definition: TGeoManager.h:600
float xmin
Definition: THbookFile.cxx:93
float ymin
Definition: THbookFile.cxx:93
float xmax
Definition: THbookFile.cxx:93
float ymax
Definition: THbookFile.cxx:93
point * points
Definition: X3DBuffer.c:22
Generic 3D primitive description class.
Definition: TBuffer3D.h:18
Int_t * fPols
Definition: TBuffer3D.h:114
UInt_t NbPnts() const
Definition: TBuffer3D.h:80
Bool_t SectionsValid(UInt_t mask) const
Definition: TBuffer3D.h:67
@ kRawSizes
Definition: TBuffer3D.h:53
void SetSectionsValid(UInt_t mask)
Definition: TBuffer3D.h:65
Int_t * fSegs
Definition: TBuffer3D.h:113
Bool_t fLocalFrame
Definition: TBuffer3D.h:90
Bool_t SetRawSizes(UInt_t reqPnts, UInt_t reqPntsCapacity, UInt_t reqSegs, UInt_t reqSegsCapacity, UInt_t reqPols, UInt_t reqPolsCapacity)
Set kRaw tessellation section of buffer with supplied sizes.
Definition: TBuffer3D.cxx:359
Double_t * fPnts
Definition: TBuffer3D.h:112
Double_t fDX
Definition: TGeoBBox.h:21
virtual Double_t DistFromOutside(const Double_t *point, const Double_t *dir, Int_t iact=1, Double_t step=TGeoShape::Big(), Double_t *safe=0) const
Compute distance from outside point to surface of the box.
Definition: TGeoBBox.cxx:430
Double_t fOrigin[3]
Definition: TGeoBBox.h:24
Double_t fDY
Definition: TGeoBBox.h:22
Double_t fDZ
Definition: TGeoBBox.h:23
virtual void FillBuffer3D(TBuffer3D &buffer, Int_t reqSections, Bool_t localFrame) const
Fills the supplied buffer, with sections in desired frame See TBuffer3D.h for explanation of sections...
Definition: TGeoBBox.cxx:1033
TGeoVolumeMulti * MakeVolumeMulti(const char *name, TGeoMedium *medium)
Make a TGeoVolumeMulti handling a list of volumes.
TObjArray * GetListOfShapes() const
Definition: TGeoManager.h:493
static Int_t ThreadId()
Translates the current thread id to an ordinal number.
Node containing an offset.
Definition: TGeoNode.h:184
Base finder class for patterns.
void SetDivIndex(Int_t index)
A polycone.
Definition: TGeoPcon.h:18
Double_t GetDphi() const
Definition: TGeoPcon.h:76
Double_t * fRmax
Definition: TGeoPcon.h:25
Double_t * fRmin
Definition: TGeoPcon.h:24
Int_t fNz
Definition: TGeoPcon.h:21
Double_t * fZ
Definition: TGeoPcon.h:26
virtual void DefineSection(Int_t snum, Double_t z, Double_t rmin, Double_t rmax)
Defines z position of a section plane, rmin and rmax at this z.
Definition: TGeoPcon.cxx:571
virtual void InspectShape() const
print shape parameters
Definition: TGeoPcon.cxx:802
Bool_t HasInsideSurface() const
Returns true when pgon has internal surface It will be only disabled when all Rmin values are 0.
Definition: TGeoPcon.cxx:1386
Double_t fPhi1
Definition: TGeoPcon.h:22
Double_t fDphi
Definition: TGeoPcon.h:23
Int_t GetNz() const
Definition: TGeoPcon.h:77
A polygone.
Definition: TGeoPgon.h:20
virtual Bool_t Contains(const Double_t *point) const
test if point is inside this shape check total z range
Definition: TGeoPgon.cxx:372
virtual void DistFromOutside_v(const Double_t *points, const Double_t *dirs, Double_t *dists, Int_t vecsize, Double_t *step) const
Compute distance from array of input points having directions specified by dirs. Store output in dist...
Definition: TGeoPgon.cxx:2154
Bool_t SliceCrossingInZ(const Double_t *point, const Double_t *dir, Int_t nphi, Int_t *iphi, Double_t *sphi, Double_t &snext, Double_t stepmax) const
Performs ray propagation between Z segments.
Definition: TGeoPgon.cxx:610
void CreateThreadData(Int_t nthreads)
Create thread data for n threads max.
Definition: TGeoPgon.cxx:107
Bool_t SliceCrossing(const Double_t *point, const Double_t *dir, Int_t nphi, Int_t *iphi, Double_t *sphi, Double_t &snext, Double_t stepmax) const
Check boundary crossing inside phi slices.
Definition: TGeoPgon.cxx:900
virtual Double_t Capacity() const
Computes capacity of the shape in [length^3].
Definition: TGeoPgon.cxx:186
virtual void SetDimensions(Double_t *param)
Set PGON dimensions starting from an array.
Definition: TGeoPgon.cxx:1938
virtual TGeoVolume * Divide(TGeoVolume *voldiv, const char *divname, Int_t iaxis, Int_t ndiv, Double_t start, Double_t step)
Divide this polygone shape belonging to volume "voldiv" into ndiv volumes called divname,...
Definition: TGeoPgon.cxx:1299
virtual Double_t DistFromInside(const Double_t *point, const Double_t *dir, Int_t iact=1, Double_t step=TGeoShape::Big(), Double_t *safe=0) const
compute distance from inside point to surface of the polygone first find out in which Z section the p...
Definition: TGeoPgon.cxx:419
virtual void ComputeNormal_v(const Double_t *points, const Double_t *dirs, Double_t *norms, Int_t vecsize)
Compute the normal for an array o points so that norm.dot.dir is positive Input: Arrays of point coor...
Definition: TGeoPgon.cxx:2137
virtual void ComputeBBox()
compute bounding box for a polygone Check if the sections are in increasing Z order
Definition: TGeoPgon.cxx:210
virtual void DistFromInside_v(const Double_t *points, const Double_t *dirs, Double_t *dists, Int_t vecsize, Double_t *step) const
Compute distance from array of input points having directions specified by dirs. Store output in dist...
Definition: TGeoPgon.cxx:2145
virtual Double_t Safety(const Double_t *point, Bool_t in=kTRUE) const
computes the closest distance from given point to this shape, according to option.
Definition: TGeoPgon.cxx:1839
virtual TBuffer3D * MakeBuffer3D() const
Creates a TBuffer3D describing this shape.
Definition: TGeoPgon.cxx:1419
Int_t fNedges
Definition: TGeoPgon.h:36
Bool_t SliceCrossingZ(const Double_t *point, const Double_t *dir, Int_t nphi, Int_t *iphi, Double_t *sphi, Double_t &snext, Double_t stepmax) const
Performs ray propagation between Z segments.
Definition: TGeoPgon.cxx:683
void LocatePhi(const Double_t *point, Int_t &ipsec) const
Locates index IPSEC of the phi sector containing POINT.
Definition: TGeoPgon.cxx:519
std::mutex fMutex
Size for the navigation data array.
Definition: TGeoPgon.h:39
TGeoPgon()
dummy ctor
Definition: TGeoPgon.cxx:125
Int_t GetNedges() const
Definition: TGeoPgon.h:82
virtual void Sizeof3D() const
fill size of this 3-D object
Definition: TGeoPgon.cxx:2085
virtual void GetBoundingCylinder(Double_t *param) const
Fill vector param[4] with the bounding cylinder parameters.
Definition: TGeoPgon.cxx:1384
virtual void GetMeshNumbers(Int_t &nvert, Int_t &nsegs, Int_t &npols) const
Returns numbers of vertices, segments and polygons composing the shape mesh.
Definition: TGeoPgon.cxx:2049
ThreadData_t & GetThreadData() const
Definition: TGeoPgon.cxx:84
virtual void SavePrimitive(std::ostream &out, Option_t *option="")
Save a primitive as a C++ statement(s) on output stream "out".
Definition: TGeoPgon.cxx:1916
virtual void SetSegsAndPols(TBuffer3D &buff) const
Fill TBuffer3D structure for segments and polygons.
Definition: TGeoPgon.cxx:1440
virtual void Contains_v(const Double_t *points, Bool_t *inside, Int_t vecsize) const
Check the inside status for each of the points in the array.
Definition: TGeoPgon.cxx:2127
Double_t SafetyToSegment(const Double_t *point, Int_t ipl, Int_t iphi, Bool_t in, Double_t safphi, Double_t safmin=TGeoShape::Big()) const
Compute safety from POINT to segment between planes ipl, ipl+1 within safmin.
Definition: TGeoPgon.cxx:1769
virtual void SetPoints(Double_t *points) const
create polygone mesh points
Definition: TGeoPgon.cxx:1963
Double_t Rpg(Double_t z, Int_t ipl, Bool_t inner, Double_t &a, Double_t &b) const
Computes projected pgon radius (inner or outer) corresponding to a given Z value.
Definition: TGeoPgon.cxx:1719
virtual Double_t DistFromOutside(const Double_t *point, const Double_t *dir, Int_t iact=1, Double_t step=TGeoShape::Big(), Double_t *safe=0) const
Compute distance from outside point to surface of the polygone.
Definition: TGeoPgon.cxx:1094
void SetSegsAndPolsNoInside(TBuffer3D &buff) const
Fill TBuffer3D structure for segments and polygons, when no inner surface exists.
Definition: TGeoPgon.cxx:1623
Bool_t SliceCrossingIn(const Double_t *point, const Double_t *dir, Int_t ipl, Int_t nphi, Int_t *iphi, Double_t *sphi, Double_t &snext, Double_t stepmax) const
Check boundary crossing inside phi slices.
Definition: TGeoPgon.cxx:764
virtual void ComputeNormal(const Double_t *point, const Double_t *dir, Double_t *norm)
Compute normal to closest surface from POINT.
Definition: TGeoPgon.cxx:277
virtual Int_t GetNmeshVertices() const
Return number of vertices of the mesh representation.
Definition: TGeoPgon.cxx:2073
std::vector< ThreadData_t * > fThreadData
Definition: TGeoPgon.h:37
Int_t fThreadSize
Navigation data per thread.
Definition: TGeoPgon.h:38
void ClearThreadData() const
Definition: TGeoPgon.cxx:92
Double_t Rproj(Double_t z, const Double_t *point, const Double_t *dir, Double_t cphi, Double_t sphi, Double_t &a, Double_t &b) const
Computes projected distance at a given Z for a given ray inside a given sector and fills coefficients...
Definition: TGeoPgon.cxx:1753
Bool_t IsCrossingSlice(const Double_t *point, const Double_t *dir, Int_t iphi, Double_t sstart, Int_t &ipl, Double_t &snext, Double_t stepmax) const
Check crossing of a given pgon slice, from a starting point inside the slice.
Definition: TGeoPgon.cxx:989
Int_t GetPhiCrossList(const Double_t *point, const Double_t *dir, Int_t istart, Double_t *sphi, Int_t *iphi, Double_t stepmax=TGeoShape::Big()) const
Mutex for thread data.
Definition: TGeoPgon.cxx:530
virtual const TBuffer3D & GetBuffer3D(Int_t reqSections, Bool_t localFrame) const
Fills a static 3D buffer and returns a reference.
Definition: TGeoPgon.cxx:2092
virtual void Safety_v(const Double_t *points, const Bool_t *inside, Double_t *safe, Int_t vecsize) const
Compute safe distance from each of the points in the input array.
Definition: TGeoPgon.cxx:2165
virtual ~TGeoPgon()
destructor
Definition: TGeoPgon.cxx:178
virtual void InspectShape() const
Inspect the PGON parameters.
Definition: TGeoPgon.cxx:1408
virtual Int_t DistancetoPrimitive(Int_t px, Int_t py)
compute closest distance from point px,py to each corner
Definition: TGeoPgon.cxx:1284
Base abstract class for all shapes.
Definition: TGeoShape.h:26
static Double_t Big()
Definition: TGeoShape.h:88
Int_t GetBasicColor() const
Get the basic color (0-7).
Definition: TGeoShape.cxx:673
void TransformPoints(Double_t *points, UInt_t NbPoints) const
Tranform a set of points (LocalToMaster)
Definition: TGeoShape.cxx:552
void SetShapeBit(UInt_t f, Bool_t set)
Equivalent of TObject::SetBit.
Definition: TGeoShape.cxx:524
static Double_t SafetyPhi(const Double_t *point, Bool_t in, Double_t phi1, Double_t phi2)
Static method to compute safety w.r.t a phi corner defined by cosines/sines of the angles phi1,...
Definition: TGeoShape.cxx:464
static Bool_t IsSameWithinTolerance(Double_t a, Double_t b)
Check if two numbers differ with less than a tolerance.
Definition: TGeoShape.cxx:326
const char * GetPointerName() const
Provide a pointer name containing uid.
Definition: TGeoShape.cxx:699
Int_t ShapeDistancetoPrimitive(Int_t numpoints, Int_t px, Int_t py) const
Returns distance to shape primitive mesh.
Definition: TGeoShape.cxx:259
virtual const char * GetName() const
Get the shape name.
Definition: TGeoShape.cxx:248
static void NormalPhi(const Double_t *point, const Double_t *dir, Double_t *norm, Double_t c1, Double_t s1, Double_t c2, Double_t s2)
Static method to compute normal to phi planes.
Definition: TGeoShape.cxx:437
static Bool_t IsCrossingSemiplane(const Double_t *point, const Double_t *dir, Double_t cphi, Double_t sphi, Double_t &snext, Double_t &rxy)
Compute distance from POINT to semiplane defined by PHI angle along DIR.
Definition: TGeoShape.cxx:297
@ kGeoClosedShape
Definition: TGeoShape.h:60
@ kGeoSavePrimitive
Definition: TGeoShape.h:65
@ kGeoPgon
Definition: TGeoShape.h:52
static Double_t Tolerance()
Definition: TGeoShape.h:91
static Bool_t IsCloseToPhi(Double_t epsil, const Double_t *point, Double_t c1, Double_t s1, Double_t c2, Double_t s2)
True if point is closer than epsil to one of the phi planes defined by c1,s1 or c2,...
Definition: TGeoShape.cxx:269
static Double_t DistFromOutsideS(const Double_t *point, const Double_t *dir, Double_t rmin, Double_t rmax, Double_t dz)
Static method to compute distance from outside point to a tube with given parameters Boundary safe al...
Definition: TGeoTube.cxx:348
Volume families.
Definition: TGeoVolume.h:252
void AddVolume(TGeoVolume *vol)
Add a volume with valid shape to the list of volumes.
TGeoVolume, TGeoVolumeMulti, TGeoVolumeAssembly are the volume classes.
Definition: TGeoVolume.h:47
void AddNodeOffset(TGeoVolume *vol, Int_t copy_no, Double_t offset=0, Option_t *option="")
Add a division node to the list of nodes.
Definition: TGeoVolume.cxx:970
TGeoMedium * GetMedium() const
Definition: TGeoVolume.h:171
void SetFinder(TGeoPatternFinder *finder)
Definition: TGeoVolume.h:229
Int_t GetNdaughters() const
Definition: TGeoVolume.h:347
TObjArray * GetNodes()
Definition: TGeoVolume.h:165
Int_t IndexOf(const TObject *obj) const
Definition: TObjArray.cxx:605
TObject * At(Int_t idx) const
Definition: TObjArray.h:166
R__ALWAYS_INLINE Bool_t TestBit(UInt_t f) const
Definition: TObject.h:187
virtual const char * ClassName() const
Returns name of class to which the object belongs.
Definition: TObject.cxx:128
void SetBit(UInt_t f, Bool_t set)
Set or unset the user status bits as specified in f.
Definition: TObject.cxx:694
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:891
virtual void Fatal(const char *method, const char *msgfmt,...) const
Issue fatal error message.
Definition: TObject.cxx:919
Basic string class.
Definition: TString.h:131
const char * Data() const
Definition: TString.h:364
TPaveText * pt
return c1
Definition: legend1.C:41
const Int_t n
Definition: legend1.C:16
return c2
Definition: legend2.C:14
double dist(Rotation3D const &r1, Rotation3D const &r2)
Definition: 3DDistances.cxx:48
Long64_t LocMin(Long64_t n, const T *a)
Return index of array with the minimum element.
Definition: TMath.h:962
Short_t Max(Short_t a, Short_t b)
Definition: TMathBase.h:212
T1 Sign(T1 a, T2 b)
Definition: TMathBase.h:165
Double_t ATan2(Double_t y, Double_t x)
Definition: TMath.h:669
Long64_t LocMax(Long64_t n, const T *a)
Return index of array with the maximum element.
Definition: TMath.h:990
constexpr Double_t E()
Base of natural log:
Definition: TMath.h:97
constexpr Double_t DegToRad()
Conversion from degree to radian:
Definition: TMath.h:82
Double_t Sqrt(Double_t x)
Definition: TMath.h:681
Short_t Min(Short_t a, Short_t b)
Definition: TMathBase.h:180
Double_t Cos(Double_t)
Definition: TMath.h:631
Double_t Sin(Double_t)
Definition: TMath.h:627
Double_t Tan(Double_t)
Definition: TMath.h:635
Long64_t BinarySearch(Long64_t n, const T *array, T value)
Definition: TMathBase.h:278
constexpr Double_t RadToDeg()
Conversion from radian to degree:
Definition: TMath.h:74
Short_t Abs(Short_t d)
Definition: TMathBase.h:120
Double_t * fDblBuffer
[fNedges+4] temporary int buffer array
Definition: TGeoPgon.h:25
ThreadData_t()
[fNedges+4] temporary double buffer array
Definition: TGeoPgon.cxx:68
~ThreadData_t()
Destructor.
Definition: TGeoPgon.cxx:76
auto * m
Definition: textangle.C:8
auto * a
Definition: textangle.C:12
#define snext(osub1, osub2)
Definition: triangle.c:1167