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