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