ROOT   Reference Guide
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/*************************************************************************
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);
38 gGeoManager->CloseGeometry();
39 gGeoManager->SetNsegments(80);
40 top->Draw();
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
75
76////////////////////////////////////////////////////////////////////////////////
77/// Destructor.
78
80{
81 delete[] fIntBuffer;
82 delete[] fDblBuffer;
83}
84
85////////////////////////////////////////////////////////////////////////////////
86
88{
91}
92
93////////////////////////////////////////////////////////////////////////////////
94
96{
97 std::lock_guard<std::mutex> guard(fMutex);
99 while (i != fThreadData.end()) {
100 delete *i;
101 ++i;
102 }
105}
106
107////////////////////////////////////////////////////////////////////////////////
109
111{
114 std::lock_guard<std::mutex> guard(fMutex);
117 for (Int_t tid = 0; tid < nthreads; tid++) {
118 if (fThreadData[tid] == nullptr) {
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;
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;
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;
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();
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;
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;
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;
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)];
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 }
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
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());
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());
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()) {
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
Definition TBuffer3D.h:67
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
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
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,...
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.
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.
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.
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
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
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
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
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