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