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