Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
TGeoCone.cxx
Go to the documentation of this file.
1// @(#)root/geom:$Id$
2// Author: Andrei Gheata 31/01/02
3// TGeoCone::Contains() and DistFromInside() 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 TGeoCone
14\ingroup Cones
15
16The cones are defined by 5 parameters:
17
18~~~{.cpp}
19TGeoCone(Double_t dz,Double_t rmin1,Double_t rmax1,
20Double_t rmin2,Double_t rmax2);
21~~~
22
23 - `rmin1:` internal radius at Z is `-dz`
24 - `rmax1:` external radius at Z is `-dz`
25 - `rmin2:` internal radius at Z is `+dz`
26 - `rmax2:` external radius at Z is `+dz`
27 - `dz:` half length in Z (a cone ranges from `-dz` to +`dz`)
28
29A cone has Z-axis as its symmetry axis.
30
31Begin_Macro
32{
33 TCanvas *c = new TCanvas("c", "c",0,0,600,600);
34 new TGeoManager("cone", "poza4");
35 TGeoMaterial *mat = new TGeoMaterial("Al", 26.98,13,2.7);
36 TGeoMedium *med = new TGeoMedium("MED",1,mat);
37 TGeoVolume *top = gGeoManager->MakeBox("TOP",med,100,100,100);
38 gGeoManager->SetTopVolume(top);
39 TGeoVolume *vol = gGeoManager->MakeCone("CONE",med, 40,10,20,35,45);
40 vol->SetLineWidth(2);
41 top->AddNode(vol,1);
42 gGeoManager->CloseGeometry();
43 gGeoManager->SetNsegments(30);
44 top->Draw();
45 TView *view = gPad->GetView();
46 if (view) view->ShowAxis();
47}
48End_Macro
49*/
50
51/** \class TGeoConeSeg
52\ingroup Cones
53
54A cone segment is a cone having a range in `phi`. The cone segment class
55derives from **`TGeoCone`**, having two extra parameters: `phi1` and
56`phi2`.
57
58~~~{.cpp}
59TGeoConeSeg(Double_t dz,Double_t rmin1,Double_t rmax1,
60Double_t rmin2,Double_t rmax2,Double_t phi1,Double_t phi2);
61~~~
62
63Parameters `phi1` and `phi2` have the same meaning and convention as for
64tube segments.
65
66Begin_Macro
67{
68 TCanvas *c = new TCanvas("c", "c",0,0,600,600);
69 new TGeoManager("coneseg", "poza5");
70 TGeoMaterial *mat = new TGeoMaterial("Al", 26.98,13,2.7);
71 TGeoMedium *med = new TGeoMedium("MED",1,mat);
72 TGeoVolume *top = gGeoManager->MakeBox("TOP",med,100,100,100);
73 gGeoManager->SetTopVolume(top);
74 TGeoVolume *vol = gGeoManager->MakeCons("CONESEG",med, 40,30,40,10,20,-30,250);
75 top->AddNode(vol,1);
76 gGeoManager->CloseGeometry();
77 gGeoManager->SetNsegments(30);
78 top->Draw();
79 TView *view = gPad->GetView();
80 if (view) view->ShowAxis();
81}
82End_Macro
83*/
84
85#include <iostream>
86
87#include "TGeoManager.h"
88#include "TGeoVolume.h"
89#include "TVirtualGeoPainter.h"
90#include "TGeoCone.h"
91#include "TBuffer3D.h"
92#include "TBuffer3DTypes.h"
93#include "TMath.h"
94
95
96////////////////////////////////////////////////////////////////////////////////
97/// Default constructor
98
100{
102 fDz = 0.0;
103 fRmin1 = 0.0;
104 fRmax1 = 0.0;
105 fRmin2 = 0.0;
106 fRmax2 = 0.0;
107}
108
109////////////////////////////////////////////////////////////////////////////////
110/// Default constructor specifying minimum and maximum radius
111
121
122////////////////////////////////////////////////////////////////////////////////
123/// Default constructor specifying minimum and maximum radius
124
126 : TGeoBBox(name, 0, 0, 0)
127{
130 if ((dz < 0) || (rmin1 < 0) || (rmax1 < 0) || (rmin2 < 0) || (rmax2 < 0)) {
132 } else
133 ComputeBBox();
134}
135
136////////////////////////////////////////////////////////////////////////////////
137/// Default constructor specifying minimum and maximum radius
138/// - param[0] = dz
139/// - param[1] = Rmin1
140/// - param[2] = Rmax1
141/// - param[3] = Rmin2
142/// - param[4] = Rmax2
143
145{
147 SetDimensions(param);
148 if ((fDz < 0) || (fRmin1 < 0) || (fRmax1 < 0) || (fRmin2 < 0) || (fRmax2 < 0))
150 else
151 ComputeBBox();
152}
153
154////////////////////////////////////////////////////////////////////////////////
155/// Computes capacity of the shape in [length^3]
156
161
162////////////////////////////////////////////////////////////////////////////////
163/// Computes capacity of the shape in [length^3]
164
166{
167 Double_t capacity = (2. * dz * TMath::Pi() / 3.) *
168 (rmax1 * rmax1 + rmax2 * rmax2 + rmax1 * rmax2 - rmin1 * rmin1 - rmin2 * rmin2 - rmin1 * rmin2);
169 return capacity;
170}
171
172////////////////////////////////////////////////////////////////////////////////
173/// destructor
174
176
177////////////////////////////////////////////////////////////////////////////////
178/// compute bounding box of the sphere
179
181{
182 TGeoBBox *box = (TGeoBBox *)this;
183 box->SetBoxDimensions(TMath::Max(fRmax1, fRmax2), TMath::Max(fRmax1, fRmax2), fDz);
184 memset(fOrigin, 0, 3 * sizeof(Double_t));
185}
186
187////////////////////////////////////////////////////////////////////////////////
188/// Compute normal to closest surface from POINT.
189
190void TGeoCone::ComputeNormal(const Double_t *point, const Double_t *dir, Double_t *norm) const
191{
192 Double_t safr, safe, phi;
193 memset(norm, 0, 3 * sizeof(Double_t));
194 phi = TMath::ATan2(point[1], point[0]);
195 Double_t cphi = TMath::Cos(phi);
196 Double_t sphi = TMath::Sin(phi);
197 Double_t ro1 = 0.5 * (fRmin1 + fRmin2);
198 Double_t tg1 = 0.5 * (fRmin2 - fRmin1) / fDz;
199 Double_t cr1 = 1. / TMath::Sqrt(1. + tg1 * tg1);
200 Double_t ro2 = 0.5 * (fRmax1 + fRmax2);
201 Double_t tg2 = 0.5 * (fRmax2 - fRmax1) / fDz;
202 Double_t cr2 = 1. / TMath::Sqrt(1. + tg2 * tg2);
203
204 Double_t r = TMath::Sqrt(point[0] * point[0] + point[1] * point[1]);
205 Double_t rin = tg1 * point[2] + ro1;
206 Double_t rout = tg2 * point[2] + ro2;
207 safe = TMath::Abs(fDz - TMath::Abs(point[2]));
208 norm[2] = 1;
209
210 safr = (ro1 > 0) ? (TMath::Abs((r - rin) * cr1)) : TGeoShape::Big();
211 if (safr < safe) {
212 safe = safr;
213 norm[0] = cr1 * cphi;
214 norm[1] = cr1 * sphi;
215 norm[2] = -tg1 * cr1;
216 }
217 safr = TMath::Abs((rout - r) * cr2);
218 if (safr < safe) {
219 norm[0] = cr2 * cphi;
220 norm[1] = cr2 * sphi;
221 norm[2] = -tg2 * cr2;
222 }
223 if (norm[0] * dir[0] + norm[1] * dir[1] + norm[2] * dir[2] < 0) {
224 norm[0] = -norm[0];
225 norm[1] = -norm[1];
226 norm[2] = -norm[2];
227 }
228}
229
230////////////////////////////////////////////////////////////////////////////////
231/// Compute normal to closest surface from POINT.
232
235{
236 Double_t safe, phi;
237 memset(norm, 0, 3 * sizeof(Double_t));
238 phi = TMath::ATan2(point[1], point[0]);
239 Double_t cphi = TMath::Cos(phi);
240 Double_t sphi = TMath::Sin(phi);
241 Double_t ro1 = 0.5 * (rmin1 + rmin2);
242 Double_t tg1 = 0.5 * (rmin2 - rmin1) / dz;
243 Double_t cr1 = 1. / TMath::Sqrt(1. + tg1 * tg1);
244 Double_t ro2 = 0.5 * (rmax1 + rmax2);
245 Double_t tg2 = 0.5 * (rmax2 - rmax1) / dz;
246 Double_t cr2 = 1. / TMath::Sqrt(1. + tg2 * tg2);
247
248 Double_t r = TMath::Sqrt(point[0] * point[0] + point[1] * point[1]);
249 Double_t rin = tg1 * point[2] + ro1;
250 Double_t rout = tg2 * point[2] + ro2;
251 safe = (ro1 > 0) ? (TMath::Abs((r - rin) * cr1)) : TGeoShape::Big();
252 norm[0] = cr1 * cphi;
253 norm[1] = cr1 * sphi;
254 norm[2] = -tg1 * cr1;
255 if (TMath::Abs((rout - r) * cr2) < safe) {
256 norm[0] = cr2 * cphi;
257 norm[1] = cr2 * sphi;
258 norm[2] = -tg2 * cr2;
259 }
260 if (norm[0] * dir[0] + norm[1] * dir[1] + norm[2] * dir[2] < 0) {
261 norm[0] = -norm[0];
262 norm[1] = -norm[1];
263 norm[2] = -norm[2];
264 }
265}
266
267////////////////////////////////////////////////////////////////////////////////
268/// test if point is inside this cone
269
271{
272 if (TMath::Abs(point[2]) > fDz)
273 return kFALSE;
274 Double_t r2 = point[0] * point[0] + point[1] * point[1];
275 Double_t rl = 0.5 * (fRmin2 * (point[2] + fDz) + fRmin1 * (fDz - point[2])) / fDz;
276 Double_t rh = 0.5 * (fRmax2 * (point[2] + fDz) + fRmax1 * (fDz - point[2])) / fDz;
277 if ((r2 < rl * rl) || (r2 > rh * rh))
278 return kFALSE;
279 return kTRUE;
280}
281
282////////////////////////////////////////////////////////////////////////////////
283/// Compute distance from inside point to surface of the cone (static)
284/// Boundary safe algorithm.
285
288{
289 if (dz <= 0)
290 return TGeoShape::Big();
291 // compute distance to surface
292 // Do Z
294 if (dir[2]) {
295 sz = (TMath::Sign(dz, dir[2]) - point[2]) / dir[2];
296 if (sz <= 0)
297 return 0.0;
298 }
299 Double_t rsq = point[0] * point[0] + point[1] * point[1];
300 Double_t zinv = 1. / dz;
301 Double_t rin = 0.5 * (rmin1 + rmin2 + (rmin2 - rmin1) * point[2] * zinv);
302 // Do Rmin
303 Double_t b, delta, zi;
304 if (rin > 0) {
305 // Protection in case point is actually outside the cone
306 if (rsq < rin * (rin + TGeoShape::Tolerance())) {
308 point[0] * dir[0] + point[1] * dir[1] + 0.5 * (rmin1 - rmin2) * dir[2] * zinv * TMath::Sqrt(rsq);
309 if (ddotn <= 0)
310 return 0.0;
311 } else {
312 TGeoCone::DistToCone(point, dir, dz, rmin1, rmin2, b, delta);
313 if (delta > 0) {
314 Double_t sr = -b - delta;
315 if (sr > 0) {
316 zi = point[2] + sr * dir[2];
317 if (TMath::Abs(zi) <= dz)
318 return TMath::Min(sz, sr);
319 }
320 sr = -b + delta;
321 if (sr > 0) {
322 zi = point[2] + sr * dir[2];
323 if (TMath::Abs(zi) <= dz)
324 return TMath::Min(sz, sr);
325 }
326 }
327 }
328 }
329 // Do Rmax
330 Double_t rout = 0.5 * (rmax1 + rmax2 + (rmax2 - rmax1) * point[2] * zinv);
331 if (rsq > rout * (rout - TGeoShape::Tolerance())) {
332 Double_t ddotn = point[0] * dir[0] + point[1] * dir[1] + 0.5 * (rmax1 - rmax2) * dir[2] * zinv * TMath::Sqrt(rsq);
333 if (ddotn >= 0)
334 return 0.0;
335 TGeoCone::DistToCone(point, dir, dz, rmax1, rmax2, b, delta);
336 if (delta < 0)
337 return 0.0;
338 Double_t sr = -b + delta;
339 if (sr < 0)
340 return sz;
341 if (TMath::Abs(-b - delta) > sr)
342 return sz;
343 zi = point[2] + sr * dir[2];
344 if (TMath::Abs(zi) <= dz)
345 return TMath::Min(sz, sr);
346 return sz;
347 }
348 TGeoCone::DistToCone(point, dir, dz, rmax1, rmax2, b, delta);
349 if (delta > 0) {
350 Double_t sr = -b - delta;
351 if (sr > 0) {
352 zi = point[2] + sr * dir[2];
353 if (TMath::Abs(zi) <= dz)
354 return TMath::Min(sz, sr);
355 }
356 sr = -b + delta;
357 if (sr > TGeoShape::Tolerance()) {
358 zi = point[2] + sr * dir[2];
359 if (TMath::Abs(zi) <= dz)
360 return TMath::Min(sz, sr);
361 }
362 }
363 return sz;
364}
365
366////////////////////////////////////////////////////////////////////////////////
367/// Compute distance from inside point to surface of the cone
368/// Boundary safe algorithm.
369
372{
373 if (iact < 3 && safe) {
374 *safe = Safety(point, kTRUE);
375 if (iact == 0)
376 return TGeoShape::Big();
377 if ((iact == 1) && (*safe > step))
378 return TGeoShape::Big();
379 }
380 // compute distance to surface
381 return TGeoCone::DistFromInsideS(point, dir, fDz, fRmin1, fRmax1, fRmin2, fRmax2);
382}
383
384////////////////////////////////////////////////////////////////////////////////
385/// Compute distance from outside point to surface of the tube
386/// Boundary safe algorithm.
387
390{
391 // compute distance to Z planes
392 if (dz <= 0)
393 return TGeoShape::Big();
395 Double_t xp, yp, zp;
396 Bool_t inz = kTRUE;
397
398 if (point[2] <= -dz) {
399 if (dir[2] <= 0)
400 return TGeoShape::Big();
401 snxt = (-dz - point[2]) / dir[2];
402 xp = point[0] + snxt * dir[0];
403 yp = point[1] + snxt * dir[1];
404 Double_t r2 = xp * xp + yp * yp;
405 if ((r2 >= rmin1 * rmin1) && (r2 <= rmax1 * rmax1))
406 return snxt;
407 inz = kFALSE;
408 } else {
409 if (point[2] >= dz) {
410 if (dir[2] >= 0)
411 return TGeoShape::Big();
412 snxt = (dz - point[2]) / dir[2];
413 xp = point[0] + snxt * dir[0];
414 yp = point[1] + snxt * dir[1];
415 Double_t r2 = xp * xp + yp * yp;
416 if ((r2 >= rmin2 * rmin2) && (r2 <= rmax2 * rmax2))
417 return snxt;
418 inz = kFALSE;
419 }
420 }
421
422 Double_t rsq = point[0] * point[0] + point[1] * point[1];
423 Double_t dzinv = 1. / dz;
424 Double_t ro1 = 0.5 * (rmin1 + rmin2);
425 Bool_t hasrmin = (ro1 > 0) ? kTRUE : kFALSE;
426 Double_t tg1 = 0.;
427 Double_t rin = 0.;
428 Bool_t inrmin = kTRUE; // r>=rmin
429 if (hasrmin) {
430 tg1 = 0.5 * (rmin2 - rmin1) * dzinv;
431 rin = ro1 + tg1 * point[2];
432 if (rin > 0 && rsq < rin * (rin - TGeoShape::Tolerance()))
433 inrmin = kFALSE;
434 }
435 Double_t ro2 = 0.5 * (rmax1 + rmax2);
436 Double_t tg2 = 0.5 * (rmax2 - rmax1) * dzinv;
437 Double_t rout = tg2 * point[2] + ro2;
439 if (rout > 0 && rsq < rout * (rout + TGeoShape::Tolerance()))
440 inrmax = kTRUE;
441 Bool_t in = inz & inrmin & inrmax;
442 Double_t b, delta;
443 // If inside cone, we are most likely on a boundary within machine precision.
444 if (in) {
446 Double_t safz = dz - TMath::Abs(point[2]); // positive
449 if (safz <= safrmin && safz <= safrmax) {
450 // on Z boundary
451 if (point[2] * dir[2] < 0)
452 return 0.0;
453 return TGeoShape::Big();
454 }
455 if (safrmax < safrmin) {
456 // on rmax boundary
457 Double_t ddotn = point[0] * dir[0] + point[1] * dir[1] - tg2 * dir[2] * r;
458 if (ddotn <= 0)
459 return 0.0;
460 return TGeoShape::Big();
461 }
462 // on rmin boundary
463 Double_t ddotn = point[0] * dir[0] + point[1] * dir[1] - tg1 * dir[2] * r;
464 if (ddotn >= 0)
465 return 0.0;
466 // we can cross (+) solution of rmin
467 TGeoCone::DistToCone(point, dir, dz, rmin1, rmin2, b, delta);
468
469 if (delta < 0)
470 return 0.0;
471 snxt = -b + delta;
472 if (snxt < 0)
473 return TGeoShape::Big();
474 if (TMath::Abs(-b - delta) > snxt)
475 return TGeoShape::Big();
476 zp = point[2] + snxt * dir[2];
477 if (TMath::Abs(zp) <= dz)
478 return snxt;
479 return TGeoShape::Big();
480 }
481
482 // compute distance to inner cone
484 if (!inrmin) {
485 // ray can cross inner cone (but not only!)
486 TGeoCone::DistToCone(point, dir, dz, rmin1, rmin2, b, delta);
487 if (delta < 0)
488 return TGeoShape::Big();
489 snxt = -b + delta;
490 if (snxt > 0) {
491 zp = point[2] + snxt * dir[2];
492 if (TMath::Abs(zp) <= dz)
493 return snxt;
494 }
495 snxt = -b - delta;
496 if (snxt > 0) {
497 zp = point[2] + snxt * dir[2];
498 if (TMath::Abs(zp) <= dz)
499 return snxt;
500 }
502 } else {
503 if (hasrmin) {
504 TGeoCone::DistToCone(point, dir, dz, rmin1, rmin2, b, delta);
505 if (delta > 0) {
506 Double_t din = -b + delta;
507 if (din > 0) {
508 zp = point[2] + din * dir[2];
509 if (TMath::Abs(zp) <= dz)
510 snxt = din;
511 }
512 }
513 }
514 }
515
516 if (inrmax)
517 return snxt;
518 // We can cross outer cone, both solutions possible
519 // compute distance to outer cone
520 TGeoCone::DistToCone(point, dir, dz, rmax1, rmax2, b, delta);
521 if (delta < 0)
522 return snxt;
523 Double_t dout = -b - delta;
524 if (dout > 0 && dout < snxt) {
525 zp = point[2] + dout * dir[2];
526 if (TMath::Abs(zp) <= dz)
527 return dout;
528 }
529 dout = -b + delta;
531 return snxt;
532 zp = point[2] + dout * dir[2];
533 if (TMath::Abs(zp) <= dz)
534 return dout;
535 return snxt;
536}
537
538////////////////////////////////////////////////////////////////////////////////
539/// compute distance from outside point to surface of the tube
540
543{
544 // compute safe radius
545 if (iact < 3 && safe) {
546 *safe = Safety(point, kFALSE);
547 if (iact == 0)
548 return TGeoShape::Big();
549 if ((iact == 1) && (*safe > step))
550 return TGeoShape::Big();
551 }
552 // Check if the bounding box is crossed within the requested distance
553 Double_t sdist = TGeoBBox::DistFromOutside(point, dir, fDX, fDY, fDZ, fOrigin, step);
554 if (sdist >= step)
555 return TGeoShape::Big();
556 // compute distance to Z planes
558}
559
560////////////////////////////////////////////////////////////////////////////////
561/// Static method to compute distance to a conical surface with :
562/// - r1, z1 - radius and Z position of lower base
563/// - r2, z2 - radius and Z position of upper base
564
566 Double_t &b, Double_t &delta)
567{
568 b = 0;
569 delta = -1.;
570 if (dz < 0)
571 return;
572 Double_t ro0 = 0.5 * (r1 + r2);
573 Double_t tz = 0.5 * (r2 - r1) / dz;
574 Double_t rsq = point[0] * point[0] + point[1] * point[1];
575 Double_t rc = ro0 + point[2] * tz;
576
577 Double_t a = dir[0] * dir[0] + dir[1] * dir[1] - tz * tz * dir[2] * dir[2];
578 b = point[0] * dir[0] + point[1] * dir[1] - tz * rc * dir[2];
579 Double_t c = rsq - rc * rc;
580
583 return;
584 b = 0.5 * c / b;
585 delta = 0.;
586 return;
587 }
588 a = 1. / a;
589 b *= a;
590 c *= a;
591 delta = b * b - c;
592 if (delta > 0) {
593 delta = TMath::Sqrt(delta);
594 } else {
595 delta = -1.;
596 }
597}
598
599////////////////////////////////////////////////////////////////////////////////
600/// compute closest distance from point px,py to each corner
601
603{
605 const Int_t numPoints = 4 * n;
606 return ShapeDistancetoPrimitive(numPoints, px, py);
607}
608
609////////////////////////////////////////////////////////////////////////////////
610/// Divide this cone shape belonging to volume "voldiv" into ndiv volumes
611/// called divname, from start position with the given step. Returns pointer
612/// to created division cell volume in case of Z divisions. For Z division
613/// creates all volumes with different shapes and returns pointer to volume that
614/// was divided. In case a wrong division axis is supplied, returns pointer to
615/// volume that was divided.
616
619{
620 TGeoShape *shape; //--- shape to be created
621 TGeoVolume *vol; //--- division volume to be created
622 TGeoVolumeMulti *vmulti; //--- generic divided volume
623 TGeoPatternFinder *finder; //--- finder to be attached
624 TString opt = ""; //--- option to be attached
625 Int_t id;
626 Double_t end = start + ndiv * step;
627 switch (iaxis) {
628 case 1: //--- R division
629 Error("Divide", "division of a cone on R not implemented");
630 return nullptr;
631 case 2: // --- Phi division
632 finder = new TGeoPatternCylPhi(voldiv, ndiv, start, end);
633 voldiv->SetFinder(finder);
634 finder->SetDivIndex(voldiv->GetNdaughters());
635 shape = new TGeoConeSeg(fDz, fRmin1, fRmax1, fRmin2, fRmax2, -step / 2, step / 2);
636 vol = new TGeoVolume(divname, shape, voldiv->GetMedium());
638 vmulti->AddVolume(vol);
639 opt = "Phi";
640 for (id = 0; id < ndiv; id++) {
641 voldiv->AddNodeOffset(vol, id, start + id * step + step / 2, opt.Data());
642 ((TGeoNodeOffset *)voldiv->GetNodes()->At(voldiv->GetNdaughters() - 1))->SetFinder(finder);
643 }
644 return vmulti;
645 case 3: //--- Z division
647 finder = new TGeoPatternZ(voldiv, ndiv, start, end);
648 voldiv->SetFinder(finder);
649 finder->SetDivIndex(voldiv->GetNdaughters());
650 for (id = 0; id < ndiv; id++) {
651 Double_t z1 = start + id * step;
652 Double_t z2 = start + (id + 1) * step;
653 Double_t rmin1n = 0.5 * (fRmin1 * (fDz - z1) + fRmin2 * (fDz + z1)) / fDz;
654 Double_t rmax1n = 0.5 * (fRmax1 * (fDz - z1) + fRmax2 * (fDz + z1)) / fDz;
655 Double_t rmin2n = 0.5 * (fRmin1 * (fDz - z2) + fRmin2 * (fDz + z2)) / fDz;
656 Double_t rmax2n = 0.5 * (fRmax1 * (fDz - z2) + fRmax2 * (fDz + z2)) / fDz;
657 shape = new TGeoCone(0.5 * step, rmin1n, rmax1n, rmin2n, rmax2n);
658 vol = new TGeoVolume(divname, shape, voldiv->GetMedium());
659 vmulti->AddVolume(vol);
660 opt = "Z";
661 voldiv->AddNodeOffset(vol, id, start + id * step + step / 2, opt.Data());
662 ((TGeoNodeOffset *)voldiv->GetNodes()->At(voldiv->GetNdaughters() - 1))->SetFinder(finder);
663 }
664 return vmulti;
665 default: Error("Divide", "Wrong axis type for division"); return nullptr;
666 }
667}
668
669////////////////////////////////////////////////////////////////////////////////
670/// Returns name of axis IAXIS.
671
673{
674 switch (iaxis) {
675 case 1: return "R";
676 case 2: return "PHI";
677 case 3: return "Z";
678 default: return "undefined";
679 }
680}
681
682////////////////////////////////////////////////////////////////////////////////
683/// Get range of shape for a given axis.
684
686{
687 xlo = 0;
688 xhi = 0;
689 Double_t dx = 0;
690 switch (iaxis) {
691 case 2:
692 xlo = 0.;
693 xhi = 360.;
694 return 360.;
695 case 3:
696 xlo = -fDz;
697 xhi = fDz;
698 dx = xhi - xlo;
699 return dx;
700 }
701 return dx;
702}
703
704////////////////////////////////////////////////////////////////////////////////
705/// Fill vector param[4] with the bounding cylinder parameters. The order
706/// is the following : Rmin, Rmax, Phi1, Phi2, dZ
707
709{
710 param[0] = TMath::Min(fRmin1, fRmin2); // Rmin
711 param[0] *= param[0];
712 param[1] = TMath::Max(fRmax1, fRmax2); // Rmax
713 param[1] *= param[1];
714 param[2] = 0.; // Phi1
715 param[3] = 360.; // Phi2
716}
717
718////////////////////////////////////////////////////////////////////////////////
719/// in case shape has some negative parameters, these has to be computed
720/// in order to fit the mother
721
723{
725 return nullptr;
726 if (!mother->TestShapeBit(kGeoCone)) {
727 Error("GetMakeRuntimeShape", "invalid mother");
728 return nullptr;
729 }
731 rmin1 = fRmin1;
732 rmax1 = fRmax1;
733 rmin2 = fRmin2;
734 rmax2 = fRmax2;
735 dz = fDz;
736 if (fDz < 0)
737 dz = ((TGeoCone *)mother)->GetDz();
738 if (fRmin1 < 0)
739 rmin1 = ((TGeoCone *)mother)->GetRmin1();
740 if (fRmax1 < 0)
741 rmax1 = ((TGeoCone *)mother)->GetRmax1();
742 if (fRmin2 < 0)
743 rmin2 = ((TGeoCone *)mother)->GetRmin2();
744 if (fRmax2 < 0)
745 rmax2 = ((TGeoCone *)mother)->GetRmax2();
746
747 return (new TGeoCone(GetName(), dz, rmin1, rmax1, rmin2, rmax2));
748}
749
750////////////////////////////////////////////////////////////////////////////////
751/// Fills array with n random points located on the line segments of the shape mesh.
752/// The output array must be provided with a length of minimum 3*npoints. Returns
753/// true if operation is implemented.
754
756{
757 if (npoints > (npoints / 2) * 2) {
758 Error("GetPointsOnSegments", "Npoints must be even number");
759 return kFALSE;
760 }
761 Bool_t hasrmin = (fRmin1 > 0 || fRmin2 > 0) ? kTRUE : kFALSE;
762 Int_t nc = 0;
763 if (hasrmin)
764 nc = (Int_t)TMath::Sqrt(0.5 * npoints);
765 else
766 nc = (Int_t)TMath::Sqrt(1. * npoints);
767 Double_t dphi = TMath::TwoPi() / nc;
768 Double_t phi = 0;
769 Int_t ntop = 0;
770 if (hasrmin)
771 ntop = npoints / 2 - nc * (nc - 1);
772 else
773 ntop = npoints - nc * (nc - 1);
774 Double_t dz = 2 * fDz / (nc - 1);
775 Double_t z = 0;
776 Int_t icrt = 0;
777 Int_t nphi = nc;
778 Double_t rmin = 0.;
779 Double_t rmax = 0.;
780 // loop z sections
781 for (Int_t i = 0; i < nc; i++) {
782 if (i == (nc - 1))
783 nphi = ntop;
784 z = -fDz + i * dz;
785 if (hasrmin)
786 rmin = 0.5 * (fRmin1 + fRmin2) + 0.5 * (fRmin2 - fRmin1) * z / fDz;
787 rmax = 0.5 * (fRmax1 + fRmax2) + 0.5 * (fRmax2 - fRmax1) * z / fDz;
788 // loop points on circle sections
789 for (Int_t j = 0; j < nphi; j++) {
790 phi = j * dphi;
791 if (hasrmin) {
792 array[icrt++] = rmin * TMath::Cos(phi);
793 array[icrt++] = rmin * TMath::Sin(phi);
794 array[icrt++] = z;
795 }
796 array[icrt++] = rmax * TMath::Cos(phi);
797 array[icrt++] = rmax * TMath::Sin(phi);
798 array[icrt++] = z;
799 }
800 }
801 return kTRUE;
802}
803
804////////////////////////////////////////////////////////////////////////////////
805/// print shape parameters
806
808{
809 printf("*** Shape %s TGeoCone ***\n", GetName());
810 printf(" dz =: %11.5f\n", fDz);
811 printf(" Rmin1 = %11.5f\n", fRmin1);
812 printf(" Rmax1 = %11.5f\n", fRmax1);
813 printf(" Rmin2 = %11.5f\n", fRmin2);
814 printf(" Rmax2 = %11.5f\n", fRmax2);
815 printf(" Bounding box:\n");
817}
818
819////////////////////////////////////////////////////////////////////////////////
820/// Creates a TBuffer3D describing *this* shape.
821/// Coordinates are in local reference frame.
822
824{
826 Int_t nbPnts = 4 * n;
827 Int_t nbSegs = 8 * n;
828 Int_t nbPols = 4 * n;
829 TBuffer3D *buff =
831
832 if (buff) {
833 SetPoints(buff->fPnts);
835 }
836 return buff;
837}
838
839////////////////////////////////////////////////////////////////////////////////
840/// Fill TBuffer3D structure for segments and polygons.
841
843{
844 Int_t i, j;
847
848 for (i = 0; i < 4; i++) {
849 for (j = 0; j < n; j++) {
850 buffer.fSegs[(i * n + j) * 3] = c;
851 buffer.fSegs[(i * n + j) * 3 + 1] = i * n + j;
852 buffer.fSegs[(i * n + j) * 3 + 2] = i * n + j + 1;
853 }
854 buffer.fSegs[(i * n + j - 1) * 3 + 2] = i * n;
855 }
856 for (i = 4; i < 6; i++) {
857 for (j = 0; j < n; j++) {
858 buffer.fSegs[(i * n + j) * 3] = c + 1;
859 buffer.fSegs[(i * n + j) * 3 + 1] = (i - 4) * n + j;
860 buffer.fSegs[(i * n + j) * 3 + 2] = (i - 2) * n + j;
861 }
862 }
863 for (i = 6; i < 8; i++) {
864 for (j = 0; j < n; j++) {
865 buffer.fSegs[(i * n + j) * 3] = c;
866 buffer.fSegs[(i * n + j) * 3 + 1] = 2 * (i - 6) * n + j;
867 buffer.fSegs[(i * n + j) * 3 + 2] = (2 * (i - 6) + 1) * n + j;
868 }
869 }
870
871 Int_t indx = 0;
872 i = 0;
873 for (j = 0; j < n; j++) {
874 indx = 6 * (i * n + j);
875 buffer.fPols[indx] = c;
876 buffer.fPols[indx + 1] = 4;
877 buffer.fPols[indx + 5] = i * n + j;
878 buffer.fPols[indx + 4] = (4 + i) * n + j;
879 buffer.fPols[indx + 3] = (2 + i) * n + j;
880 buffer.fPols[indx + 2] = (4 + i) * n + j + 1;
881 }
882 buffer.fPols[indx + 2] = (4 + i) * n;
883 i = 1;
884 for (j = 0; j < n; j++) {
885 indx = 6 * (i * n + j);
886 buffer.fPols[indx] = c;
887 buffer.fPols[indx + 1] = 4;
888 buffer.fPols[indx + 2] = i * n + j;
889 buffer.fPols[indx + 3] = (4 + i) * n + j;
890 buffer.fPols[indx + 4] = (2 + i) * n + j;
891 buffer.fPols[indx + 5] = (4 + i) * n + j + 1;
892 }
893 buffer.fPols[indx + 5] = (4 + i) * n;
894 i = 2;
895 for (j = 0; j < n; j++) {
896 indx = 6 * (i * n + j);
897 buffer.fPols[indx] = c + i;
898 buffer.fPols[indx + 1] = 4;
899 buffer.fPols[indx + 2] = (i - 2) * 2 * n + j;
900 buffer.fPols[indx + 3] = (4 + i) * n + j;
901 buffer.fPols[indx + 4] = ((i - 2) * 2 + 1) * n + j;
902 buffer.fPols[indx + 5] = (4 + i) * n + j + 1;
903 }
904 buffer.fPols[indx + 5] = (4 + i) * n;
905 i = 3;
906 for (j = 0; j < n; j++) {
907 indx = 6 * (i * n + j);
908 buffer.fPols[indx] = c + i;
909 buffer.fPols[indx + 1] = 4;
910 buffer.fPols[indx + 5] = (i - 2) * 2 * n + j;
911 buffer.fPols[indx + 4] = (4 + i) * n + j;
912 buffer.fPols[indx + 3] = ((i - 2) * 2 + 1) * n + j;
913 buffer.fPols[indx + 2] = (4 + i) * n + j + 1;
914 }
915 buffer.fPols[indx + 2] = (4 + i) * n;
916}
917
918////////////////////////////////////////////////////////////////////////////////
919/// computes the closest distance from given point to this shape, according
920/// to option. The matching point on the shape is stored in spoint.
921
923{
924 Double_t saf[4];
925 Double_t r = TMath::Sqrt(point[0] * point[0] + point[1] * point[1]);
926 saf[0] = TGeoShape::SafetySeg(r, point[2], fRmin1, -fDz, fRmax1, -fDz, !in);
927 saf[1] = TGeoShape::SafetySeg(r, point[2], fRmax2, fDz, fRmin2, fDz, !in);
928 saf[2] = TGeoShape::SafetySeg(r, point[2], fRmin2, fDz, fRmin1, -fDz, !in);
929 saf[3] = TGeoShape::SafetySeg(r, point[2], fRmax1, -fDz, fRmax2, fDz, !in);
931 if (safety > 1.E20)
932 safety = 0.;
933 return safety;
934}
935
936////////////////////////////////////////////////////////////////////////////////
937/// computes the closest distance from given point to this shape, according
938/// to option. The matching point on the shape is stored in spoint.
939
942{
943 Double_t saf[4];
944 Double_t r = TMath::Sqrt(point[0] * point[0] + point[1] * point[1]);
945 // Double_t rin = tg1*point[2]+ro1;
946 // Double_t rout = tg2*point[2]+ro2;
947 switch (skipz) {
948 case 1: // skip lower Z plane
949 saf[0] = TGeoShape::Big();
950 saf[1] = TGeoShape::SafetySeg(r, point[2], rmax2, dz, rmin2, dz, !in);
951 break;
952 case 2: // skip upper Z plane
953 saf[0] = TGeoShape::SafetySeg(r, point[2], rmin1, -dz, rmax1, -dz, !in);
954 saf[1] = TGeoShape::Big();
955 break;
956 case 3: // skip both
957 saf[0] = saf[1] = TGeoShape::Big();
958 break;
959 default:
960 saf[0] = TGeoShape::SafetySeg(r, point[2], rmin1, -dz, rmax1, -dz, !in);
961 saf[1] = TGeoShape::SafetySeg(r, point[2], rmax2, dz, rmin2, dz, !in);
962 }
963 // Safety to inner part
964 if (rmin1 > 0 || rmin2 > 0)
965 saf[2] = TGeoShape::SafetySeg(r, point[2], rmin2, dz, rmin1, -dz, !in);
966 else
967 saf[2] = TGeoShape::Big();
968 saf[3] = TGeoShape::SafetySeg(r, point[2], rmax1, -dz, rmax2, dz, !in);
969 return saf[TMath::LocMin(4, saf)];
970}
971
972////////////////////////////////////////////////////////////////////////////////
973/// Save a primitive as a C++ statement(s) on output stream "out".
974
975void TGeoCone::SavePrimitive(std::ostream &out, Option_t * /*option*/ /*= ""*/)
976{
978 return;
979 out << " // Shape: " << GetName() << " type: " << ClassName() << std::endl;
980 out << " dz = " << fDz << ";" << std::endl;
981 out << " rmin1 = " << fRmin1 << ";" << std::endl;
982 out << " rmax1 = " << fRmax1 << ";" << std::endl;
983 out << " rmin2 = " << fRmin2 << ";" << std::endl;
984 out << " rmax2 = " << fRmax2 << ";" << std::endl;
985 out << " TGeoShape *" << GetPointerName() << " = new TGeoCone(\"" << GetName()
986 << "\", dz,rmin1,rmax1,rmin2,rmax2);" << std::endl;
988}
989
990////////////////////////////////////////////////////////////////////////////////
991/// Set cone dimensions.
992
994{
995 if (rmin1 >= 0) {
996 if (rmax1 > 0) {
997 if (rmin1 <= rmax1) {
998 // normal rmin/rmax
999 fRmin1 = rmin1;
1000 fRmax1 = rmax1;
1001 } else {
1002 fRmin1 = rmax1;
1003 fRmax1 = rmin1;
1004 Warning("SetConeDimensions", "rmin1>rmax1 Switch rmin1<->rmax1");
1006 }
1007 } else {
1008 // run-time
1009 fRmin1 = rmin1;
1010 fRmax1 = rmax1;
1011 }
1012 } else {
1013 // run-time
1014 fRmin1 = rmin1;
1015 fRmax1 = rmax1;
1016 }
1017 if (rmin2 >= 0) {
1018 if (rmax2 > 0) {
1019 if (rmin2 <= rmax2) {
1020 // normal rmin/rmax
1021 fRmin2 = rmin2;
1022 fRmax2 = rmax2;
1023 } else {
1024 fRmin2 = rmax2;
1025 fRmax2 = rmin2;
1026 Warning("SetConeDimensions", "rmin2>rmax2 Switch rmin2<->rmax2");
1028 }
1029 } else {
1030 // run-time
1031 fRmin2 = rmin2;
1032 fRmax2 = rmax2;
1033 }
1034 } else {
1035 // run-time
1036 fRmin2 = rmin2;
1037 fRmax2 = rmax2;
1038 }
1039
1040 fDz = dz;
1041}
1042
1043////////////////////////////////////////////////////////////////////////////////
1044/// Set cone dimensions from an array.
1045
1047{
1048 Double_t dz = param[0];
1049 Double_t rmin1 = param[1];
1050 Double_t rmax1 = param[2];
1051 Double_t rmin2 = param[3];
1052 Double_t rmax2 = param[4];
1054}
1055
1056////////////////////////////////////////////////////////////////////////////////
1057/// Create cone mesh points.
1058
1060{
1061 Double_t dz, phi, dphi;
1062 Int_t j, n;
1063
1065 dphi = 360. / n;
1066 dz = fDz;
1067 Int_t indx = 0;
1068
1069 if (points) {
1070 for (j = 0; j < n; j++) {
1071 phi = j * dphi * TMath::DegToRad();
1072 points[indx++] = fRmin1 * TMath::Cos(phi);
1073 points[indx++] = fRmin1 * TMath::Sin(phi);
1074 points[indx++] = -dz;
1075 }
1076
1077 for (j = 0; j < n; j++) {
1078 phi = j * dphi * TMath::DegToRad();
1079 points[indx++] = fRmax1 * TMath::Cos(phi);
1080 points[indx++] = fRmax1 * TMath::Sin(phi);
1081 points[indx++] = -dz;
1082 }
1083
1084 for (j = 0; j < n; j++) {
1085 phi = j * dphi * TMath::DegToRad();
1086 points[indx++] = fRmin2 * TMath::Cos(phi);
1087 points[indx++] = fRmin2 * TMath::Sin(phi);
1088 points[indx++] = dz;
1089 }
1090
1091 for (j = 0; j < n; j++) {
1092 phi = j * dphi * TMath::DegToRad();
1093 points[indx++] = fRmax2 * TMath::Cos(phi);
1094 points[indx++] = fRmax2 * TMath::Sin(phi);
1095 points[indx++] = dz;
1096 }
1097 }
1098}
1099
1100////////////////////////////////////////////////////////////////////////////////
1101/// Create cone mesh points.
1102
1104{
1105 Double_t dz, phi, dphi;
1106 Int_t j, n;
1107
1109 dphi = 360. / n;
1110 dz = fDz;
1111 Int_t indx = 0;
1112
1113 if (points) {
1114 for (j = 0; j < n; j++) {
1115 phi = j * dphi * TMath::DegToRad();
1116 points[indx++] = fRmin1 * TMath::Cos(phi);
1117 points[indx++] = fRmin1 * TMath::Sin(phi);
1118 points[indx++] = -dz;
1119 }
1120
1121 for (j = 0; j < n; j++) {
1122 phi = j * dphi * TMath::DegToRad();
1123 points[indx++] = fRmax1 * TMath::Cos(phi);
1124 points[indx++] = fRmax1 * TMath::Sin(phi);
1125 points[indx++] = -dz;
1126 }
1127
1128 for (j = 0; j < n; j++) {
1129 phi = j * dphi * TMath::DegToRad();
1130 points[indx++] = fRmin2 * TMath::Cos(phi);
1131 points[indx++] = fRmin2 * TMath::Sin(phi);
1132 points[indx++] = dz;
1133 }
1134
1135 for (j = 0; j < n; j++) {
1136 phi = j * dphi * TMath::DegToRad();
1137 points[indx++] = fRmax2 * TMath::Cos(phi);
1138 points[indx++] = fRmax2 * TMath::Sin(phi);
1139 points[indx++] = dz;
1140 }
1141 }
1142}
1143
1144////////////////////////////////////////////////////////////////////////////////
1145/// Returns numbers of vertices, segments and polygons composing the shape mesh.
1146
1148{
1150 nvert = n * 4;
1151 nsegs = n * 8;
1152 npols = n * 4;
1153}
1154
1155////////////////////////////////////////////////////////////////////////////////
1156/// Return number of vertices of the mesh representation
1157
1159{
1161 Int_t numPoints = n * 4;
1162 return numPoints;
1163}
1164
1165////////////////////////////////////////////////////////////////////////////////
1166/// Fill size of this 3-D object
1167
1168void TGeoCone::Sizeof3D() const {}
1169
1170////////////////////////////////////////////////////////////////////////////////
1171/// Fills a static 3D buffer and returns a reference.
1172
1174{
1175 static TBuffer3D buffer(TBuffer3DTypes::kGeneric);
1176
1178
1181 Int_t nbPnts = 4 * n;
1182 Int_t nbSegs = 8 * n;
1183 Int_t nbPols = 4 * n;
1184 if (buffer.SetRawSizes(nbPnts, 3 * nbPnts, nbSegs, 3 * nbSegs, nbPols, 6 * nbPols)) {
1186 }
1187 }
1188
1189 // TODO: Can we push this as common down to TGeoShape?
1191 SetPoints(buffer.fPnts);
1192 if (!buffer.fLocalFrame) {
1193 TransformPoints(buffer.fPnts, buffer.NbPnts());
1194 }
1195
1196 SetSegsAndPols(buffer);
1198 }
1199
1200 return buffer;
1201}
1202
1203////////////////////////////////////////////////////////////////////////////////
1204/// Check the inside status for each of the points in the array.
1205/// Input: Array of point coordinates + vector size
1206/// Output: Array of Booleans for the inside of each point
1207
1209{
1210 for (Int_t i = 0; i < vecsize; i++)
1211 inside[i] = Contains(&points[3 * i]);
1212}
1213
1214////////////////////////////////////////////////////////////////////////////////
1215/// Compute the normal for an array o points so that norm.dot.dir is positive
1216/// Input: Arrays of point coordinates and directions + vector size
1217/// Output: Array of normal directions
1218
1220{
1221 for (Int_t i = 0; i < vecsize; i++)
1222 ComputeNormal(&points[3 * i], &dirs[3 * i], &norms[3 * i]);
1223}
1224
1225////////////////////////////////////////////////////////////////////////////////
1226/// Compute distance from array of input points having directions specified by dirs. Store output in dists
1227
1229 Double_t *step) const
1230{
1231 for (Int_t i = 0; i < vecsize; i++)
1232 dists[i] = DistFromInside(&points[3 * i], &dirs[3 * i], 3, step[i]);
1233}
1234
1235////////////////////////////////////////////////////////////////////////////////
1236/// Compute distance from array of input points having directions specified by dirs. Store output in dists
1237
1239 Double_t *step) const
1240{
1241 for (Int_t i = 0; i < vecsize; i++)
1242 dists[i] = DistFromOutside(&points[3 * i], &dirs[3 * i], 3, step[i]);
1243}
1244
1245////////////////////////////////////////////////////////////////////////////////
1246/// Compute safe distance from each of the points in the input array.
1247/// Input: Array of point coordinates, array of statuses for these points, size of the arrays
1248/// Output: Safety values
1249
1251{
1252 for (Int_t i = 0; i < vecsize; i++)
1253 safe[i] = Safety(&points[3 * i], inside[i]);
1254}
1255
1256
1257////////////////////////////////////////////////////////////////////////////////
1258/// Default constructor
1259
1261 : TGeoCone(), fPhi1(0.), fPhi2(0.), fS1(0.), fC1(0.), fS2(0.), fC2(0.), fSm(0.), fCm(0.), fCdfi(0.)
1262{
1264 fPhi1 = fPhi2 = 0.0;
1265}
1266
1267////////////////////////////////////////////////////////////////////////////////
1268/// Default constructor specifying minimum and maximum radius
1269
1271 Double_t phi2)
1273 fPhi1(0.),
1274 fPhi2(0.),
1275 fS1(0.),
1276 fC1(0.),
1277 fS2(0.),
1278 fC2(0.),
1279 fSm(0.),
1280 fCm(0.),
1281 fCdfi(0.)
1282
1283{
1286 ComputeBBox();
1287}
1288
1289////////////////////////////////////////////////////////////////////////////////
1290/// Default constructor specifying minimum and maximum radius
1291
1295 fPhi1(0.),
1296 fPhi2(0.),
1297 fS1(0.),
1298 fC1(0.),
1299 fS2(0.),
1300 fC2(0.),
1301 fSm(0.),
1302 fCm(0.),
1303 fCdfi(0.)
1304{
1307 ComputeBBox();
1308}
1309
1310////////////////////////////////////////////////////////////////////////////////
1311/// Default constructor specifying minimum and maximum radius
1312/// - param[0] = dz
1313/// - param[1] = Rmin1
1314/// - param[2] = Rmax1
1315/// - param[3] = Rmin2
1316/// - param[4] = Rmax2
1317/// - param[5] = phi1
1318/// - param[6] = phi2
1319
1321 : TGeoCone(0, 0, 0, 0, 0), fPhi1(0.), fPhi2(0.), fS1(0.), fC1(0.), fS2(0.), fC2(0.), fSm(0.), fCm(0.), fCdfi(0.)
1322{
1324 SetDimensions(param);
1325 ComputeBBox();
1326}
1327
1328////////////////////////////////////////////////////////////////////////////////
1329/// destructor
1330
1332
1333////////////////////////////////////////////////////////////////////////////////
1334/// Function called after streaming an object of this class.
1335
1340
1341////////////////////////////////////////////////////////////////////////////////
1342/// Init frequently used trigonometric values
1343
1345{
1348 fC1 = TMath::Cos(phi1);
1349 fS1 = TMath::Sin(phi1);
1350 fC2 = TMath::Cos(phi2);
1351 fS2 = TMath::Sin(phi2);
1352 Double_t fio = 0.5 * (phi1 + phi2);
1353 fCm = TMath::Cos(fio);
1354 fSm = TMath::Sin(fio);
1355 Double_t dfi = 0.5 * (phi2 - phi1);
1356 fCdfi = TMath::Cos(dfi);
1357}
1358
1359////////////////////////////////////////////////////////////////////////////////
1360/// Computes capacity of the shape in [length^3]
1361
1366
1367////////////////////////////////////////////////////////////////////////////////
1368/// Computes capacity of the shape in [length^3]
1369
1377
1378////////////////////////////////////////////////////////////////////////////////
1379/// compute bounding box of the tube segment
1380
1382{
1386
1387 Double_t xc[4];
1388 Double_t yc[4];
1389 xc[0] = rmax * fC1;
1390 yc[0] = rmax * fS1;
1391 xc[1] = rmax * fC2;
1392 yc[1] = rmax * fS2;
1393 xc[2] = rmin * fC1;
1394 yc[2] = rmin * fS1;
1395 xc[3] = rmin * fC2;
1396 yc[3] = rmin * fS2;
1397
1398 Double_t xmin = xc[TMath::LocMin(4, &xc[0])];
1399 Double_t xmax = xc[TMath::LocMax(4, &xc[0])];
1400 Double_t ymin = yc[TMath::LocMin(4, &yc[0])];
1401 Double_t ymax = yc[TMath::LocMax(4, &yc[0])];
1402
1403 Double_t dp = fPhi2 - fPhi1;
1404 Double_t ddp = -fPhi1;
1405 if (ddp < 0)
1406 ddp += 360;
1407 if (ddp <= dp)
1408 xmax = rmax;
1409 ddp = 90 - fPhi1;
1410 if (ddp < 0)
1411 ddp += 360;
1412 if (ddp <= dp)
1413 ymax = rmax;
1414 ddp = 180 - fPhi1;
1415 if (ddp < 0)
1416 ddp += 360;
1417 if (ddp <= dp)
1418 xmin = -rmax;
1419 ddp = 270 - fPhi1;
1420 if (ddp < 0)
1421 ddp += 360;
1422 if (ddp <= dp)
1423 ymin = -rmax;
1424 fOrigin[0] = (xmax + xmin) / 2;
1425 fOrigin[1] = (ymax + ymin) / 2;
1426 fOrigin[2] = 0;
1427 fDX = (xmax - xmin) / 2;
1428 fDY = (ymax - ymin) / 2;
1429 fDZ = fDz;
1430}
1431
1432////////////////////////////////////////////////////////////////////////////////
1433/// Compute normal to closest surface from POINT.
1434
1435void TGeoConeSeg::ComputeNormal(const Double_t *point, const Double_t *dir, Double_t *norm) const
1436{
1437 Double_t saf[3];
1438 Double_t ro1 = 0.5 * (fRmin1 + fRmin2);
1439 Double_t tg1 = 0.5 * (fRmin2 - fRmin1) / fDz;
1440 Double_t cr1 = 1. / TMath::Sqrt(1. + tg1 * tg1);
1441 Double_t ro2 = 0.5 * (fRmax1 + fRmax2);
1442 Double_t tg2 = 0.5 * (fRmax2 - fRmax1) / fDz;
1443 Double_t cr2 = 1. / TMath::Sqrt(1. + tg2 * tg2);
1444
1445 Double_t r = TMath::Sqrt(point[0] * point[0] + point[1] * point[1]);
1446 Double_t rin = tg1 * point[2] + ro1;
1447 Double_t rout = tg2 * point[2] + ro2;
1448 saf[0] = TMath::Abs(fDz - TMath::Abs(point[2]));
1449 saf[1] = (ro1 > 0) ? (TMath::Abs((r - rin) * cr1)) : TGeoShape::Big();
1450 saf[2] = TMath::Abs((rout - r) * cr2);
1451 Int_t i = TMath::LocMin(3, saf);
1452 if (((fPhi2 - fPhi1) < 360.) && TGeoShape::IsCloseToPhi(saf[i], point, fC1, fS1, fC2, fS2)) {
1453 TGeoShape::NormalPhi(point, dir, norm, fC1, fS1, fC2, fS2);
1454 return;
1455 }
1456 if (i == 0) {
1457 norm[0] = norm[1] = 0.;
1458 norm[2] = TMath::Sign(1., dir[2]);
1459 return;
1460 }
1461
1462 Double_t phi = TMath::ATan2(point[1], point[0]);
1463 Double_t cphi = TMath::Cos(phi);
1464 Double_t sphi = TMath::Sin(phi);
1465
1466 if (i == 1) {
1467 norm[0] = cr1 * cphi;
1468 norm[1] = cr1 * sphi;
1469 norm[2] = -tg1 * cr1;
1470 } else {
1471 norm[0] = cr2 * cphi;
1472 norm[1] = cr2 * sphi;
1473 norm[2] = -tg2 * cr2;
1474 }
1475
1476 if (norm[0] * dir[0] + norm[1] * dir[1] + norm[2] * dir[2] < 0) {
1477 norm[0] = -norm[0];
1478 norm[1] = -norm[1];
1479 norm[2] = -norm[2];
1480 }
1481}
1482
1483////////////////////////////////////////////////////////////////////////////////
1484/// Compute normal to closest surface from POINT.
1485
1489{
1490 Double_t saf[2];
1491 Double_t ro1 = 0.5 * (rmin1 + rmin2);
1492 Double_t tg1 = 0.5 * (rmin2 - rmin1) / dz;
1493 Double_t cr1 = 1. / TMath::Sqrt(1. + tg1 * tg1);
1494 Double_t ro2 = 0.5 * (rmax1 + rmax2);
1495 Double_t tg2 = 0.5 * (rmax2 - rmax1) / dz;
1496 Double_t cr2 = 1. / TMath::Sqrt(1. + tg2 * tg2);
1497
1498 Double_t r = TMath::Sqrt(point[0] * point[0] + point[1] * point[1]);
1499 Double_t rin = tg1 * point[2] + ro1;
1500 Double_t rout = tg2 * point[2] + ro2;
1501 saf[0] = (ro1 > 0) ? (TMath::Abs((r - rin) * cr1)) : TGeoShape::Big();
1502 saf[1] = TMath::Abs((rout - r) * cr2);
1503 Int_t i = TMath::LocMin(2, saf);
1504 if (TGeoShape::IsCloseToPhi(saf[i], point, c1, s1, c2, s2)) {
1505 TGeoShape::NormalPhi(point, dir, norm, c1, s1, c2, s2);
1506 return;
1507 }
1508
1509 Double_t phi = TMath::ATan2(point[1], point[0]);
1510 Double_t cphi = TMath::Cos(phi);
1511 Double_t sphi = TMath::Sin(phi);
1512
1513 if (i == 0) {
1514 norm[0] = cr1 * cphi;
1515 norm[1] = cr1 * sphi;
1516 norm[2] = -tg1 * cr1;
1517 } else {
1518 norm[0] = cr2 * cphi;
1519 norm[1] = cr2 * sphi;
1520 norm[2] = -tg2 * cr2;
1521 }
1522
1523 if (norm[0] * dir[0] + norm[1] * dir[1] + norm[2] * dir[2] < 0) {
1524 norm[0] = -norm[0];
1525 norm[1] = -norm[1];
1526 norm[2] = -norm[2];
1527 }
1528}
1529
1530////////////////////////////////////////////////////////////////////////////////
1531/// test if point is inside this sphere
1532
1534{
1535 if (!TGeoCone::Contains(point))
1536 return kFALSE;
1538 if (dphi >= 360.)
1539 return kTRUE;
1540 Double_t phi = TMath::ATan2(point[1], point[0]) * TMath::RadToDeg();
1541 if (phi < 0)
1542 phi += 360.;
1543 Double_t ddp = phi - fPhi1;
1544 if (ddp < 0)
1545 ddp += 360.;
1546 // if (ddp > 360) ddp-=360;
1547 if (ddp > dphi)
1548 return kFALSE;
1549 return kTRUE;
1550}
1551
1552////////////////////////////////////////////////////////////////////////////////
1553/// Static method to compute distance to a conical surface with :
1554/// - r1, z1 - radius and Z position of lower base
1555/// - r2, z2 - radius and Z position of upper base
1556/// - phi1, phi2 - phi limits
1557
1560{
1561 Double_t dz = z2 - z1;
1562 if (dz <= 0) {
1563 return TGeoShape::Big();
1564 }
1565
1566 Double_t dphi = phi2 - phi1;
1568 if (dphi >= 360.)
1569 hasphi = kFALSE;
1570 if (dphi < 0)
1571 dphi += 360.;
1572 // printf("phi1=%f phi2=%f dphi=%f\n", phi1, phi2, dphi);
1573
1574 Double_t ro0 = 0.5 * (r1 + r2);
1575 Double_t fz = (r2 - r1) / dz;
1576 Double_t r0sq = point[0] * point[0] + point[1] * point[1];
1577 Double_t rc = ro0 + fz * (point[2] - 0.5 * (z1 + z2));
1578
1579 Double_t a = dir[0] * dir[0] + dir[1] * dir[1] - fz * fz * dir[2] * dir[2];
1580 Double_t b = point[0] * dir[0] + point[1] * dir[1] - fz * rc * dir[2];
1581 Double_t c = r0sq - rc * rc;
1582
1583 if (a == 0)
1584 return TGeoShape::Big();
1585 a = 1. / a;
1586 b *= a;
1587 c *= a;
1588 Double_t delta = b * b - c;
1589 if (delta < 0)
1590 return TGeoShape::Big();
1591 delta = TMath::Sqrt(delta);
1592
1593 Double_t snxt = -b - delta;
1594 Double_t ptnew[3];
1595 Double_t ddp, phi;
1596 if (snxt > 0) {
1597 // check Z range
1598 ptnew[2] = point[2] + snxt * dir[2];
1599 if (((ptnew[2] - z1) * (ptnew[2] - z2)) < 0) {
1600 // check phi range
1601 if (!hasphi)
1602 return snxt;
1603 ptnew[0] = point[0] + snxt * dir[0];
1604 ptnew[1] = point[1] + snxt * dir[1];
1605 phi = TMath::ATan2(ptnew[1], ptnew[0]) * TMath::RadToDeg();
1606 if (phi < 0)
1607 phi += 360.;
1608 ddp = phi - phi1;
1609 if (ddp < 0)
1610 ddp += 360.;
1611 // printf("snxt1=%f phi=%f ddp=%f\n", snxt, phi, ddp);
1612 if (ddp <= dphi)
1613 return snxt;
1614 }
1615 }
1616 snxt = -b + delta;
1617 if (snxt > 0) {
1618 // check Z range
1619 ptnew[2] = point[2] + snxt * dir[2];
1620 if (((ptnew[2] - z1) * (ptnew[2] - z2)) < 0) {
1621 // check phi range
1622 if (!hasphi)
1623 return snxt;
1624 ptnew[0] = point[0] + snxt * dir[0];
1625 ptnew[1] = point[1] + snxt * dir[1];
1626 phi = TMath::ATan2(ptnew[1], ptnew[0]) * TMath::RadToDeg();
1627 if (phi < 0)
1628 phi += 360.;
1629 ddp = phi - phi1;
1630 if (ddp < 0)
1631 ddp += 360.;
1632 // printf("snxt2=%f phi=%f ddp=%f\n", snxt, phi, ddp);
1633 if (ddp <= dphi)
1634 return snxt;
1635 }
1636 }
1637 return TGeoShape::Big();
1638}
1639
1640////////////////////////////////////////////////////////////////////////////////
1641/// compute distance from inside point to surface of the tube segment
1642
1646{
1647 if (dz <= 0)
1648 return TGeoShape::Big();
1649 // Do Z
1651 if (scone <= 0)
1652 return 0.0;
1654 Double_t rsq = point[0] * point[0] + point[1] * point[1];
1656 Double_t cpsi = point[0] * cm + point[1] * sm;
1657 if (cpsi > r * cdfi + TGeoShape::Tolerance()) {
1658 sfmin = TGeoShape::DistToPhiMin(point, dir, s1, c1, s2, c2, sm, cm);
1659 return TMath::Min(scone, sfmin);
1660 }
1661 // Point on the phi boundary or outside
1662 // which one: phi1 or phi2
1663 Double_t ddotn, xi, yi;
1664 if (TMath::Abs(point[1] - s1 * r) < TMath::Abs(point[1] - s2 * r)) {
1665 ddotn = s1 * dir[0] - c1 * dir[1];
1666 if (ddotn >= 0)
1667 return 0.0;
1668 ddotn = -s2 * dir[0] + c2 * dir[1];
1669 if (ddotn <= 0)
1670 return scone;
1671 sfmin = s2 * point[0] - c2 * point[1];
1672 if (sfmin <= 0)
1673 return scone;
1674 sfmin /= ddotn;
1675 if (sfmin >= scone)
1676 return scone;
1677 xi = point[0] + sfmin * dir[0];
1678 yi = point[1] + sfmin * dir[1];
1679 if (yi * cm - xi * sm < 0)
1680 return scone;
1681 return sfmin;
1682 }
1683 ddotn = -s2 * dir[0] + c2 * dir[1];
1684 if (ddotn >= 0)
1685 return 0.0;
1686 ddotn = s1 * dir[0] - c1 * dir[1];
1687 if (ddotn <= 0)
1688 return scone;
1689 sfmin = -s1 * point[0] + c1 * point[1];
1690 if (sfmin <= 0)
1691 return scone;
1692 sfmin /= ddotn;
1693 if (sfmin >= scone)
1694 return scone;
1695 xi = point[0] + sfmin * dir[0];
1696 yi = point[1] + sfmin * dir[1];
1697 if (yi * cm - xi * sm > 0)
1698 return scone;
1699 return sfmin;
1700}
1701
1702////////////////////////////////////////////////////////////////////////////////
1703/// compute distance from inside point to surface of the tube segment
1704
1707{
1708 if (iact < 3 && safe) {
1710 if (iact == 0)
1711 return TGeoShape::Big();
1712 if ((iact == 1) && (*safe > step))
1713 return TGeoShape::Big();
1714 }
1715 if ((fPhi2 - fPhi1) >= 360.)
1716 return TGeoCone::DistFromInsideS(point, dir, fDz, fRmin1, fRmax1, fRmin2, fRmax2);
1717
1718 // compute distance to surface
1720 fCdfi);
1721}
1722
1723////////////////////////////////////////////////////////////////////////////////
1724/// compute distance from outside point to surface of arbitrary tube
1725
1729{
1730 if (dz <= 0)
1731 return TGeoShape::Big();
1732 Double_t r2, cpsi;
1733 // check Z planes
1734 Double_t xi, yi, zi;
1735 Double_t b, delta;
1736 zi = dz - TMath::Abs(point[2]);
1737 Double_t rin, rout;
1740 Bool_t in = kFALSE;
1741 Bool_t inz = (zi < 0) ? kFALSE : kTRUE;
1742 if (!inz) {
1743 if (point[2] * dir[2] >= 0)
1744 return TGeoShape::Big();
1745 s = -zi / TMath::Abs(dir[2]);
1746 xi = point[0] + s * dir[0];
1747 yi = point[1] + s * dir[1];
1748 r2 = xi * xi + yi * yi;
1749 if (dir[2] > 0) {
1750 rin = rmin1;
1751 rout = rmax1;
1752 } else {
1753 rin = rmin2;
1754 rout = rmax2;
1755 }
1756 if ((rin * rin <= r2) && (r2 <= rout * rout)) {
1757 cpsi = xi * cm + yi * sm;
1758 if (cpsi >= (cdfi * TMath::Sqrt(r2)))
1759 return s;
1760 }
1761 }
1762 Double_t zinv = 1. / dz;
1763 Double_t rsq = point[0] * point[0] + point[1] * point[1];
1765 Double_t ro1 = 0.5 * (rmin1 + rmin2);
1766 Bool_t hasrmin = (ro1 > 0) ? kTRUE : kFALSE;
1767 Double_t tg1 = 0.0;
1769 rin = 0.0;
1770 if (hasrmin) {
1771 tg1 = 0.5 * (rmin2 - rmin1) * zinv;
1772 rin = ro1 + tg1 * point[2];
1773 if (rsq > rin * (rin - TGeoShape::Tolerance()))
1774 inrmin = kTRUE;
1775 } else {
1776 inrmin = kTRUE;
1777 }
1778 Double_t ro2 = 0.5 * (rmax1 + rmax2);
1779 Double_t tg2 = 0.5 * (rmax2 - rmax1) * zinv;
1780 rout = ro2 + tg2 * point[2];
1782 if (r < rout + TGeoShape::Tolerance())
1783 inrmax = kTRUE;
1785 cpsi = point[0] * cm + point[1] * sm;
1786 if (cpsi > r * cdfi - TGeoShape::Tolerance())
1787 inphi = kTRUE;
1788 in = inz & inrmin & inrmax & inphi;
1789 // If inside, we are most likely on a boundary within machine precision.
1790 if (in) {
1791 Double_t safphi = (cpsi - r * cdfi) * TMath::Sqrt(1. - cdfi * cdfi);
1794 // check if on Z boundaries
1795 if (zi < safrmax && zi < safrmin && zi < safphi) {
1796 if (point[2] * dir[2] < 0)
1797 return 0.0;
1798 return TGeoShape::Big();
1799 }
1800 // check if on Rmax boundary
1801 if (safrmax < safrmin && safrmax < safphi) {
1802 Double_t ddotn = point[0] * dir[0] + point[1] * dir[1] - tg2 * dir[2] * r;
1803 if (ddotn <= 0)
1804 return 0.0;
1805 return TGeoShape::Big();
1806 }
1807 // check if on phi boundary
1808 if (safphi < safrmin) {
1809 // We may cross again a phi of rmin boundary
1810 // check first if we are on phi1 or phi2
1811 Double_t un;
1812 if (point[0] * c1 + point[1] * s1 > point[0] * c2 + point[1] * s2) {
1813 un = dir[0] * s1 - dir[1] * c1;
1814 if (un < 0)
1815 return 0.0;
1816 if (cdfi >= 0)
1817 return TGeoShape::Big();
1818 un = -dir[0] * s2 + dir[1] * c2;
1819 if (un < 0) {
1820 s = -point[0] * s2 + point[1] * c2;
1821 if (s > 0) {
1822 s /= (-un);
1823 zi = point[2] + s * dir[2];
1824 if (TMath::Abs(zi) <= dz) {
1825 xi = point[0] + s * dir[0];
1826 yi = point[1] + s * dir[1];
1827 if ((yi * cm - xi * sm) > 0) {
1828 r2 = xi * xi + yi * yi;
1829 rin = ro1 + tg1 * zi;
1830 rout = ro2 + tg2 * zi;
1831 if ((rin * rin <= r2) && (rout * rout >= r2))
1832 return s;
1833 }
1834 }
1835 }
1836 }
1837 } else {
1838 un = -dir[0] * s2 + dir[1] * c2;
1839 if (un < 0)
1840 return 0.0;
1841 if (cdfi >= 0)
1842 return TGeoShape::Big();
1843 un = dir[0] * s1 - dir[1] * c1;
1844 if (un < 0) {
1845 s = point[0] * s1 - point[1] * c1;
1846 if (s > 0) {
1847 s /= (-un);
1848 zi = point[2] + s * dir[2];
1849 if (TMath::Abs(zi) <= dz) {
1850 xi = point[0] + s * dir[0];
1851 yi = point[1] + s * dir[1];
1852 if ((yi * cm - xi * sm) < 0) {
1853 r2 = xi * xi + yi * yi;
1854 rin = ro1 + tg1 * zi;
1855 rout = ro2 + tg2 * zi;
1856 if ((rin * rin <= r2) && (rout * rout >= r2))
1857 return s;
1858 }
1859 }
1860 }
1861 }
1862 }
1863 // We may also cross rmin, second solution coming from outside
1864 Double_t ddotn = point[0] * dir[0] + point[1] * dir[1] - tg1 * dir[2] * r;
1865 if (ddotn >= 0)
1866 return TGeoShape::Big();
1867 if (cdfi >= 0)
1868 return TGeoShape::Big();
1869 TGeoCone::DistToCone(point, dir, dz, rmin1, rmin2, b, delta);
1870 if (delta < 0)
1871 return TGeoShape::Big();
1872 snxt = -b - delta;
1873 if (snxt < 0)
1874 return TGeoShape::Big();
1875 snxt = -b + delta;
1876 zi = point[2] + snxt * dir[2];
1877 if (TMath::Abs(zi) > dz)
1878 return TGeoShape::Big();
1879 xi = point[0] + snxt * dir[0];
1880 yi = point[1] + snxt * dir[1];
1881 r2 = xi * xi + yi * yi;
1882 cpsi = xi * cm + yi * sm;
1883 if (cpsi >= (cdfi * TMath::Sqrt(r2)))
1884 return snxt;
1885 return TGeoShape::Big();
1886 }
1887 // We are on rmin boundary: we may cross again rmin or a phi facette
1888 Double_t ddotn = point[0] * dir[0] + point[1] * dir[1] - tg1 * dir[2] * r;
1889 if (ddotn >= 0)
1890 return 0.0;
1891 TGeoCone::DistToCone(point, dir, dz, rmin1, rmin2, b, delta);
1892 if (delta < 0)
1893 return 0.0;
1894 snxt = -b + delta;
1895 if (snxt < 0)
1896 return TGeoShape::Big();
1897 if (TMath::Abs(-b - delta) > snxt)
1898 return TGeoShape::Big();
1899 zi = point[2] + snxt * dir[2];
1900 if (TMath::Abs(zi) > dz)
1901 return TGeoShape::Big();
1902 // OK, we cross rmin at snxt - check if within phi range
1903 xi = point[0] + snxt * dir[0];
1904 yi = point[1] + snxt * dir[1];
1905 r2 = xi * xi + yi * yi;
1906 cpsi = xi * cm + yi * sm;
1907 if (cpsi >= (cdfi * TMath::Sqrt(r2)))
1908 return snxt;
1909 // we cross rmin in the phi gap - we may cross a phi facette
1910 if (cdfi >= 0)
1911 return TGeoShape::Big();
1912 Double_t un = -dir[0] * s1 + dir[1] * c1;
1913 if (un > 0) {
1914 s = point[0] * s1 - point[1] * c1;
1915 if (s >= 0) {
1916 s /= un;
1917 zi = point[2] + s * dir[2];
1918 if (TMath::Abs(zi) <= dz) {
1919 xi = point[0] + s * dir[0];
1920 yi = point[1] + s * dir[1];
1921 if ((yi * cm - xi * sm) <= 0) {
1922 r2 = xi * xi + yi * yi;
1923 rin = ro1 + tg1 * zi;
1924 rout = ro2 + tg2 * zi;
1925 if ((rin * rin <= r2) && (rout * rout >= r2))
1926 return s;
1927 }
1928 }
1929 }
1930 }
1931 un = dir[0] * s2 - dir[1] * c2;
1932 if (un > 0) {
1933 s = (point[1] * c2 - point[0] * s2) / un;
1934 if (s >= 0) {
1935 zi = point[2] + s * dir[2];
1936 if (TMath::Abs(zi) <= dz) {
1937 xi = point[0] + s * dir[0];
1938 yi = point[1] + s * dir[1];
1939 if ((yi * cm - xi * sm) >= 0) {
1940 r2 = xi * xi + yi * yi;
1941 rin = ro1 + tg1 * zi;
1942 rout = ro2 + tg2 * zi;
1943 if ((rin * rin <= r2) && (rout * rout >= r2))
1944 return s;
1945 }
1946 }
1947 }
1948 }
1949 return TGeoShape::Big();
1950 }
1951
1952 // The point is really outside
1954 if (!inrmax) {
1955 // check crossing with outer cone
1956 TGeoCone::DistToCone(point, dir, dz, rmax1, rmax2, b, delta);
1957 if (delta >= 0) {
1958 s = -b - delta;
1959 if (s > 0) {
1960 zi = point[2] + s * dir[2];
1961 if (TMath::Abs(zi) <= dz) {
1962 xi = point[0] + s * dir[0];
1963 yi = point[1] + s * dir[1];
1964 r2 = xi * xi + yi * yi;
1965 cpsi = xi * cm + yi * sm;
1966 if (cpsi >= (cdfi * TMath::Sqrt(r2)))
1967 return s; // rmax crossing
1968 }
1969 }
1970 s = -b + delta;
1971 if (s > 0) {
1972 zi = point[2] + s * dir[2];
1973 if (TMath::Abs(zi) <= dz) {
1974 xi = point[0] + s * dir[0];
1975 yi = point[1] + s * dir[1];
1976 r2 = xi * xi + yi * yi;
1977 cpsi = xi * cm + yi * sm;
1978 if (cpsi >= (cdfi * TMath::Sqrt(r2)))
1979 sr1 = s;
1980 }
1981 }
1982 }
1983 }
1984 // check crossing with inner cone
1986 TGeoCone::DistToCone(point, dir, dz, rmin1, rmin2, b, delta);
1987 if (delta >= 0) {
1988 s = -b - delta;
1989 if (s > 0) {
1990 zi = point[2] + s * dir[2];
1991 if (TMath::Abs(zi) <= dz) {
1992 xi = point[0] + s * dir[0];
1993 yi = point[1] + s * dir[1];
1994 r2 = xi * xi + yi * yi;
1995 cpsi = xi * cm + yi * sm;
1996 if (cpsi >= (cdfi * TMath::Sqrt(r2)))
1997 sr2 = s;
1998 }
1999 }
2000 if (sr2 > 1E10) {
2001 s = -b + delta;
2002 if (s > 0) {
2003 zi = point[2] + s * dir[2];
2004 if (TMath::Abs(zi) <= dz) {
2005 xi = point[0] + s * dir[0];
2006 yi = point[1] + s * dir[1];
2007 r2 = xi * xi + yi * yi;
2008 cpsi = xi * cm + yi * sm;
2009 if (cpsi >= (cdfi * TMath::Sqrt(r2)))
2010 sr2 = s;
2011 }
2012 }
2013 }
2014 }
2015 snxt = TMath::Min(sr1, sr2);
2016 // Check phi crossing
2017 s = TGeoShape::DistToPhiMin(point, dir, s1, c1, s2, c2, sm, cm, kFALSE);
2018 if (s > snxt)
2019 return snxt;
2020 zi = point[2] + s * dir[2];
2021 if (TMath::Abs(zi) > dz)
2022 return snxt;
2023 xi = point[0] + s * dir[0];
2024 yi = point[1] + s * dir[1];
2025 r2 = xi * xi + yi * yi;
2026 rout = ro2 + tg2 * zi;
2027 if (r2 > rout * rout)
2028 return snxt;
2029 rin = ro1 + tg1 * zi;
2030 if (r2 >= rin * rin)
2031 return s; // phi crossing
2032 return snxt;
2033}
2034
2035////////////////////////////////////////////////////////////////////////////////
2036/// compute distance from outside point to surface of the tube
2037
2039 Double_t *safe) const
2040{
2041 // compute safe radius
2042 if (iact < 3 && safe) {
2043 *safe = Safety(point, kFALSE);
2044 if (iact == 0)
2045 return TGeoShape::Big();
2046 if ((iact == 1) && (*safe > step))
2047 return TGeoShape::Big();
2048 }
2049 // Check if the bounding box is crossed within the requested distance
2050 Double_t sdist = TGeoBBox::DistFromOutside(point, dir, fDX, fDY, fDZ, fOrigin, step);
2051 if (sdist >= step)
2052 return TGeoShape::Big();
2053 if ((fPhi2 - fPhi1) >= 360.)
2054 return TGeoCone::DistFromOutsideS(point, dir, fDz, fRmin1, fRmax1, fRmin2, fRmax2);
2056 fCdfi);
2057}
2058
2059////////////////////////////////////////////////////////////////////////////////
2060/// compute closest distance from point px,py to each corner
2061
2063{
2065 const Int_t numPoints = 4 * n;
2066 return ShapeDistancetoPrimitive(numPoints, px, py);
2067}
2068
2069////////////////////////////////////////////////////////////////////////////////
2070/// Divide this cone segment shape belonging to volume "voldiv" into ndiv volumes
2071/// called divname, from start position with the given step. Returns pointer
2072/// to created division cell volume in case of Z divisions. For Z division
2073/// creates all volumes with different shapes and returns pointer to volume that
2074/// was divided. In case a wrong division axis is supplied, returns pointer to
2075/// volume that was divided.
2076
2077TGeoVolume *
2079{
2080 TGeoShape *shape; //--- shape to be created
2081 TGeoVolume *vol; //--- division volume to be created
2082 TGeoVolumeMulti *vmulti; //--- generic divided volume
2083 TGeoPatternFinder *finder; //--- finder to be attached
2084 TString opt = ""; //--- option to be attached
2085 Double_t dphi;
2086 Int_t id;
2087 Double_t end = start + ndiv * step;
2088 switch (iaxis) {
2089 case 1: //--- R division
2090 Error("Divide", "division of a cone segment on R not implemented");
2091 return nullptr;
2092 case 2: //--- Phi division
2093 dphi = fPhi2 - fPhi1;
2094 if (dphi < 0)
2095 dphi += 360.;
2096 finder = new TGeoPatternCylPhi(voldiv, ndiv, start, end);
2097 voldiv->SetFinder(finder);
2098 finder->SetDivIndex(voldiv->GetNdaughters());
2099 shape = new TGeoConeSeg(fDz, fRmin1, fRmax1, fRmin2, fRmax2, -step / 2, step / 2);
2100 vol = new TGeoVolume(divname, shape, voldiv->GetMedium());
2102 vmulti->AddVolume(vol);
2103 opt = "Phi";
2104 for (id = 0; id < ndiv; id++) {
2105 voldiv->AddNodeOffset(vol, id, start + id * step + step / 2, opt.Data());
2106 ((TGeoNodeOffset *)voldiv->GetNodes()->At(voldiv->GetNdaughters() - 1))->SetFinder(finder);
2107 }
2108 return vmulti;
2109 case 3: //--- Z division
2110 finder = new TGeoPatternZ(voldiv, ndiv, start, end);
2112 voldiv->SetFinder(finder);
2113 finder->SetDivIndex(voldiv->GetNdaughters());
2114 for (id = 0; id < ndiv; id++) {
2115 Double_t z1 = start + id * step;
2116 Double_t z2 = start + (id + 1) * step;
2117 Double_t rmin1n = 0.5 * (fRmin1 * (fDz - z1) + fRmin2 * (fDz + z1)) / fDz;
2118 Double_t rmax1n = 0.5 * (fRmax1 * (fDz - z1) + fRmax2 * (fDz + z1)) / fDz;
2119 Double_t rmin2n = 0.5 * (fRmin1 * (fDz - z2) + fRmin2 * (fDz + z2)) / fDz;
2120 Double_t rmax2n = 0.5 * (fRmax1 * (fDz - z2) + fRmax2 * (fDz + z2)) / fDz;
2121 shape = new TGeoConeSeg(step / 2, rmin1n, rmax1n, rmin2n, rmax2n, fPhi1, fPhi2);
2122 vol = new TGeoVolume(divname, shape, voldiv->GetMedium());
2123 vmulti->AddVolume(vol);
2124 opt = "Z";
2125 voldiv->AddNodeOffset(vol, id, start + id * step + step / 2, opt.Data());
2126 ((TGeoNodeOffset *)voldiv->GetNodes()->At(voldiv->GetNdaughters() - 1))->SetFinder(finder);
2127 }
2128 return vmulti;
2129 default: Error("Divide", "Wrong axis type for division"); return nullptr;
2130 }
2131}
2132
2133////////////////////////////////////////////////////////////////////////////////
2134/// Get range of shape for a given axis.
2135
2137{
2138 xlo = 0;
2139 xhi = 0;
2140 Double_t dx = 0;
2141 switch (iaxis) {
2142 case 2:
2143 xlo = fPhi1;
2144 xhi = fPhi2;
2145 dx = xhi - xlo;
2146 return dx;
2147 case 3:
2148 xlo = -fDz;
2149 xhi = fDz;
2150 dx = xhi - xlo;
2151 return dx;
2152 }
2153 return dx;
2154}
2155
2156////////////////////////////////////////////////////////////////////////////////
2157/// Fill vector param[4] with the bounding cylinder parameters. The order
2158/// is the following : Rmin, Rmax, Phi1, Phi2
2159
2161{
2162 param[0] = TMath::Min(fRmin1, fRmin2); // Rmin
2163 param[0] *= param[0];
2164 param[1] = TMath::Max(fRmax1, fRmax2); // Rmax
2165 param[1] *= param[1];
2166 param[2] = (fPhi1 < 0) ? (fPhi1 + 360.) : fPhi1; // Phi1
2167 param[3] = fPhi2; // Phi2
2168 while (param[3] < param[2])
2169 param[3] += 360.;
2170}
2171
2172////////////////////////////////////////////////////////////////////////////////
2173/// in case shape has some negative parameters, these has to be computed
2174/// in order to fit the mother
2175
2177{
2179 return nullptr;
2180 if (!mother->TestShapeBit(kGeoConeSeg)) {
2181 Error("GetMakeRuntimeShape", "invalid mother");
2182 return nullptr;
2183 }
2185 rmin1 = fRmin1;
2186 rmax1 = fRmax1;
2187 rmin2 = fRmin2;
2188 rmax2 = fRmax2;
2189 dz = fDz;
2190 if (fDz < 0)
2191 dz = ((TGeoCone *)mother)->GetDz();
2192 if (fRmin1 < 0)
2193 rmin1 = ((TGeoCone *)mother)->GetRmin1();
2194 if ((fRmax1 < 0) || (fRmax1 < fRmin1))
2195 rmax1 = ((TGeoCone *)mother)->GetRmax1();
2196 if (fRmin2 < 0)
2197 rmin2 = ((TGeoCone *)mother)->GetRmin2();
2198 if ((fRmax2 < 0) || (fRmax2 < fRmin2))
2199 rmax2 = ((TGeoCone *)mother)->GetRmax2();
2200
2201 return (new TGeoConeSeg(GetName(), dz, rmin1, rmax1, rmin2, rmax2, fPhi1, fPhi2));
2202}
2203
2204////////////////////////////////////////////////////////////////////////////////
2205/// print shape parameters
2206
2208{
2209 printf("*** Shape %s: TGeoConeSeg ***\n", GetName());
2210 printf(" dz = %11.5f\n", fDz);
2211 printf(" Rmin1 = %11.5f\n", fRmin1);
2212 printf(" Rmax1 = %11.5f\n", fRmax1);
2213 printf(" Rmin2 = %11.5f\n", fRmin2);
2214 printf(" Rmax2 = %11.5f\n", fRmax2);
2215 printf(" phi1 = %11.5f\n", fPhi1);
2216 printf(" phi2 = %11.5f\n", fPhi2);
2217 printf(" Bounding box:\n");
2219}
2220
2221///////////////////////////////////////////////////////////////////////////////
2222/// Creates a TBuffer3D describing *this* shape.
2223/// Coordinates are in local reference frame.
2224
2226{
2228 Int_t nbPnts = 4 * n;
2229 Int_t nbSegs = 2 * nbPnts;
2230 Int_t nbPols = nbPnts - 2;
2231 TBuffer3D *buff =
2233
2234 if (buff) {
2235 SetPoints(buff->fPnts);
2237 }
2238
2239 return buff;
2240}
2241
2242////////////////////////////////////////////////////////////////////////////////
2243/// Fill TBuffer3D structure for segments and polygons.
2244
2246{
2247 Int_t i, j;
2249 Int_t c = GetBasicColor();
2250
2251 memset(buffer.fSegs, 0, buffer.NbSegs() * 3 * sizeof(Int_t));
2252 for (i = 0; i < 4; i++) {
2253 for (j = 1; j < n; j++) {
2254 buffer.fSegs[(i * n + j - 1) * 3] = c;
2255 buffer.fSegs[(i * n + j - 1) * 3 + 1] = i * n + j - 1;
2256 buffer.fSegs[(i * n + j - 1) * 3 + 2] = i * n + j;
2257 }
2258 }
2259 for (i = 4; i < 6; i++) {
2260 for (j = 0; j < n; j++) {
2261 buffer.fSegs[(i * n + j) * 3] = c + 1;
2262 buffer.fSegs[(i * n + j) * 3 + 1] = (i - 4) * n + j;
2263 buffer.fSegs[(i * n + j) * 3 + 2] = (i - 2) * n + j;
2264 }
2265 }
2266 for (i = 6; i < 8; i++) {
2267 for (j = 0; j < n; j++) {
2268 buffer.fSegs[(i * n + j) * 3] = c;
2269 buffer.fSegs[(i * n + j) * 3 + 1] = 2 * (i - 6) * n + j;
2270 buffer.fSegs[(i * n + j) * 3 + 2] = (2 * (i - 6) + 1) * n + j;
2271 }
2272 }
2273
2274 Int_t indx = 0;
2275 memset(buffer.fPols, 0, buffer.NbPols() * 6 * sizeof(Int_t));
2276 i = 0;
2277 for (j = 0; j < n - 1; j++) {
2278 buffer.fPols[indx++] = c;
2279 buffer.fPols[indx++] = 4;
2280 buffer.fPols[indx++] = (4 + i) * n + j + 1;
2281 buffer.fPols[indx++] = (2 + i) * n + j;
2282 buffer.fPols[indx++] = (4 + i) * n + j;
2283 buffer.fPols[indx++] = i * n + j;
2284 }
2285 i = 1;
2286 for (j = 0; j < n - 1; j++) {
2287 buffer.fPols[indx++] = c;
2288 buffer.fPols[indx++] = 4;
2289 buffer.fPols[indx++] = i * n + j;
2290 buffer.fPols[indx++] = (4 + i) * n + j;
2291 buffer.fPols[indx++] = (2 + i) * n + j;
2292 buffer.fPols[indx++] = (4 + i) * n + j + 1;
2293 }
2294 i = 2;
2295 for (j = 0; j < n - 1; j++) {
2296 buffer.fPols[indx++] = c + i;
2297 buffer.fPols[indx++] = 4;
2298 buffer.fPols[indx++] = (i - 2) * 2 * n + j;
2299 buffer.fPols[indx++] = (4 + i) * n + j;
2300 buffer.fPols[indx++] = ((i - 2) * 2 + 1) * n + j;
2301 buffer.fPols[indx++] = (4 + i) * n + j + 1;
2302 }
2303 i = 3;
2304 for (j = 0; j < n - 1; j++) {
2305 buffer.fPols[indx++] = c + i;
2306 buffer.fPols[indx++] = 4;
2307 buffer.fPols[indx++] = (4 + i) * n + j + 1;
2308 buffer.fPols[indx++] = ((i - 2) * 2 + 1) * n + j;
2309 buffer.fPols[indx++] = (4 + i) * n + j;
2310 buffer.fPols[indx++] = (i - 2) * 2 * n + j;
2311 }
2312 buffer.fPols[indx++] = c + 2;
2313 buffer.fPols[indx++] = 4;
2314 buffer.fPols[indx++] = 6 * n;
2315 buffer.fPols[indx++] = 4 * n;
2316 buffer.fPols[indx++] = 7 * n;
2317 buffer.fPols[indx++] = 5 * n;
2318 buffer.fPols[indx++] = c + 2;
2319 buffer.fPols[indx++] = 4;
2320 buffer.fPols[indx++] = 6 * n - 1;
2321 buffer.fPols[indx++] = 8 * n - 1;
2322 buffer.fPols[indx++] = 5 * n - 1;
2323 buffer.fPols[indx++] = 7 * n - 1;
2324}
2325
2326////////////////////////////////////////////////////////////////////////////////
2327/// computes the closest distance from given point to this shape, according
2328/// to option. The matching point on the shape is stored in spoint.
2329
2331{
2332 Double_t safe = TGeoCone::Safety(point, in);
2333 if ((fPhi2 - fPhi1) >= 360.)
2334 return safe;
2336 if (in)
2337 return TMath::Min(safe, safphi);
2338 if (safe > 1.E10)
2339 return safphi;
2340 return TMath::Max(safe, safphi);
2341}
2342
2343////////////////////////////////////////////////////////////////////////////////
2344/// Static method to compute the closest distance from given point to this shape.
2345
2348{
2350 if ((phi2 - phi1) >= 360.)
2351 return safe;
2353 if (in)
2354 return TMath::Min(safe, safphi);
2355 if (safe > 1.E10)
2356 return safphi;
2357 return TMath::Max(safe, safphi);
2358}
2359
2360////////////////////////////////////////////////////////////////////////////////
2361/// Save a primitive as a C++ statement(s) on output stream "out".
2362
2363void TGeoConeSeg::SavePrimitive(std::ostream &out, Option_t * /*option*/ /*= ""*/)
2364{
2366 return;
2367 out << " // Shape: " << GetName() << " type: " << ClassName() << std::endl;
2368 out << " dz = " << fDz << ";" << std::endl;
2369 out << " rmin1 = " << fRmin1 << ";" << std::endl;
2370 out << " rmax1 = " << fRmax1 << ";" << std::endl;
2371 out << " rmin2 = " << fRmin2 << ";" << std::endl;
2372 out << " rmax2 = " << fRmax2 << ";" << std::endl;
2373 out << " phi1 = " << fPhi1 << ";" << std::endl;
2374 out << " phi2 = " << fPhi2 << ";" << std::endl;
2375 out << " TGeoShape *" << GetPointerName() << " = new TGeoConeSeg(\"" << GetName()
2376 << "\", dz,rmin1,rmax1,rmin2,rmax2,phi1,phi2);" << std::endl;
2378}
2379
2380////////////////////////////////////////////////////////////////////////////////
2381/// Set dimensions of the cone segment.
2382
2385{
2386 fDz = dz;
2387 fRmin1 = rmin1;
2388 fRmax1 = rmax1;
2389 fRmin2 = rmin2;
2390 fRmax2 = rmax2;
2391 fPhi1 = phi1;
2392 while (fPhi1 < 0)
2393 fPhi1 += 360.;
2394 fPhi2 = phi2;
2395 while (fPhi2 <= fPhi1)
2396 fPhi2 += 360.;
2398 Fatal("SetConsDimensions", "In shape %s invalid phi1=%g, phi2=%g\n", GetName(), fPhi1, fPhi2);
2400}
2401
2402////////////////////////////////////////////////////////////////////////////////
2403/// Set dimensions of the cone segment from an array.
2404
2406{
2407 Double_t dz = param[0];
2408 Double_t rmin1 = param[1];
2409 Double_t rmax1 = param[2];
2410 Double_t rmin2 = param[3];
2411 Double_t rmax2 = param[4];
2412 Double_t phi1 = param[5];
2413 Double_t phi2 = param[6];
2415}
2416
2417////////////////////////////////////////////////////////////////////////////////
2418/// Create cone segment mesh points.
2419
2421{
2422 Int_t j, n;
2423 Float_t dphi, phi, phi1, phi2, dz;
2424
2425 n = gGeoManager->GetNsegments() + 1;
2426 dz = fDz;
2427 phi1 = fPhi1;
2428 phi2 = fPhi2;
2429
2430 dphi = (phi2 - phi1) / (n - 1);
2431
2432 Int_t indx = 0;
2433
2434 if (points) {
2435 for (j = 0; j < n; j++) {
2436 phi = (fPhi1 + j * dphi) * TMath::DegToRad();
2437 points[indx++] = fRmin1 * TMath::Cos(phi);
2438 points[indx++] = fRmin1 * TMath::Sin(phi);
2439 points[indx++] = -dz;
2440 }
2441 for (j = 0; j < n; j++) {
2442 phi = (fPhi1 + j * dphi) * TMath::DegToRad();
2443 points[indx++] = fRmax1 * TMath::Cos(phi);
2444 points[indx++] = fRmax1 * TMath::Sin(phi);
2445 points[indx++] = -dz;
2446 }
2447 for (j = 0; j < n; j++) {
2448 phi = (fPhi1 + j * dphi) * TMath::DegToRad();
2449 points[indx++] = fRmin2 * TMath::Cos(phi);
2450 points[indx++] = fRmin2 * TMath::Sin(phi);
2451 points[indx++] = dz;
2452 }
2453 for (j = 0; j < n; j++) {
2454 phi = (fPhi1 + j * dphi) * TMath::DegToRad();
2455 points[indx++] = fRmax2 * TMath::Cos(phi);
2456 points[indx++] = fRmax2 * TMath::Sin(phi);
2457 points[indx++] = dz;
2458 }
2459 }
2460}
2461
2462////////////////////////////////////////////////////////////////////////////////
2463/// Create cone segment mesh points.
2464
2466{
2467 Int_t j, n;
2468 Float_t dphi, phi, phi1, phi2, dz;
2469
2470 n = gGeoManager->GetNsegments() + 1;
2471 dz = fDz;
2472 phi1 = fPhi1;
2473 phi2 = fPhi2;
2474
2475 dphi = (phi2 - phi1) / (n - 1);
2476
2477 Int_t indx = 0;
2478
2479 if (points) {
2480 for (j = 0; j < n; j++) {
2481 phi = (fPhi1 + j * dphi) * TMath::DegToRad();
2482 points[indx++] = fRmin1 * TMath::Cos(phi);
2483 points[indx++] = fRmin1 * TMath::Sin(phi);
2484 points[indx++] = -dz;
2485 }
2486 for (j = 0; j < n; j++) {
2487 phi = (fPhi1 + j * dphi) * TMath::DegToRad();
2488 points[indx++] = fRmax1 * TMath::Cos(phi);
2489 points[indx++] = fRmax1 * TMath::Sin(phi);
2490 points[indx++] = -dz;
2491 }
2492 for (j = 0; j < n; j++) {
2493 phi = (fPhi1 + j * dphi) * TMath::DegToRad();
2494 points[indx++] = fRmin2 * TMath::Cos(phi);
2495 points[indx++] = fRmin2 * TMath::Sin(phi);
2496 points[indx++] = dz;
2497 }
2498 for (j = 0; j < n; j++) {
2499 phi = (fPhi1 + j * dphi) * TMath::DegToRad();
2500 points[indx++] = fRmax2 * TMath::Cos(phi);
2501 points[indx++] = fRmax2 * TMath::Sin(phi);
2502 points[indx++] = dz;
2503 }
2504 }
2505}
2506
2507////////////////////////////////////////////////////////////////////////////////
2508/// Returns numbers of vertices, segments and polygons composing the shape mesh.
2509
2511{
2513 nvert = n * 4;
2514 nsegs = n * 8;
2515 npols = n * 4 - 2;
2516}
2517
2518////////////////////////////////////////////////////////////////////////////////
2519/// Return number of vertices of the mesh representation
2520
2522{
2524 Int_t numPoints = n * 4;
2525 return numPoints;
2526}
2527
2528////////////////////////////////////////////////////////////////////////////////
2529/// Fill size of this 3-D object
2530
2532
2533////////////////////////////////////////////////////////////////////////////////
2534/// Fills a static 3D buffer and returns a reference.
2535
2537{
2538 static TBuffer3D buffer(TBuffer3DTypes::kGeneric);
2539
2541
2544 Int_t nbPnts = 4 * n;
2545 Int_t nbSegs = 2 * nbPnts;
2546 Int_t nbPols = nbPnts - 2;
2547 if (buffer.SetRawSizes(nbPnts, 3 * nbPnts, nbSegs, 3 * nbSegs, nbPols, 6 * nbPols)) {
2549 }
2550 }
2552 SetPoints(buffer.fPnts);
2553 if (!buffer.fLocalFrame) {
2554 TransformPoints(buffer.fPnts, buffer.NbPnts());
2555 }
2556
2557 SetSegsAndPols(buffer);
2559 }
2560
2561 return buffer;
2562}
2563
2564////////////////////////////////////////////////////////////////////////////////
2565/// Fills array with n random points located on the line segments of the shape mesh.
2566/// The output array must be provided with a length of minimum 3*npoints. Returns
2567/// true if operation is implemented.
2568
2570{
2571 if (npoints > (npoints / 2) * 2) {
2572 Error("GetPointsOnSegments", "Npoints must be even number");
2573 return kFALSE;
2574 }
2575 Int_t nc = (Int_t)TMath::Sqrt(0.5 * npoints);
2576 Double_t dphi = (fPhi2 - fPhi1) * TMath::DegToRad() / (nc - 1);
2577 Double_t phi = 0;
2579 Int_t ntop = npoints / 2 - nc * (nc - 1);
2580 Double_t dz = 2 * fDz / (nc - 1);
2581 Double_t z = 0;
2582 Double_t rmin = 0.;
2583 Double_t rmax = 0.;
2584 Int_t icrt = 0;
2585 Int_t nphi = nc;
2586 // loop z sections
2587 for (Int_t i = 0; i < nc; i++) {
2588 if (i == (nc - 1)) {
2589 nphi = ntop;
2590 dphi = (fPhi2 - fPhi1) * TMath::DegToRad() / (nphi - 1);
2591 }
2592 z = -fDz + i * dz;
2593 rmin = 0.5 * (fRmin1 + fRmin2) + 0.5 * (fRmin2 - fRmin1) * z / fDz;
2594 rmax = 0.5 * (fRmax1 + fRmax2) + 0.5 * (fRmax2 - fRmax1) * z / fDz;
2595 // loop points on circle sections
2596 for (Int_t j = 0; j < nphi; j++) {
2597 phi = phi1 + j * dphi;
2598 array[icrt++] = rmin * TMath::Cos(phi);
2599 array[icrt++] = rmin * TMath::Sin(phi);
2600 array[icrt++] = z;
2601 array[icrt++] = rmax * TMath::Cos(phi);
2602 array[icrt++] = rmax * TMath::Sin(phi);
2603 array[icrt++] = z;
2604 }
2605 }
2606 return kTRUE;
2607}
2608
2609////////////////////////////////////////////////////////////////////////////////
2610/// Check the inside status for each of the points in the array.
2611/// Input: Array of point coordinates + vector size
2612/// Output: Array of Booleans for the inside of each point
2613
2615{
2616 for (Int_t i = 0; i < vecsize; i++)
2617 inside[i] = Contains(&points[3 * i]);
2618}
2619
2620////////////////////////////////////////////////////////////////////////////////
2621/// Compute the normal for an array o points so that norm.dot.dir is positive
2622/// Input: Arrays of point coordinates and directions + vector size
2623/// Output: Array of normal directions
2624
2626{
2627 for (Int_t i = 0; i < vecsize; i++)
2628 ComputeNormal(&points[3 * i], &dirs[3 * i], &norms[3 * i]);
2629}
2630
2631////////////////////////////////////////////////////////////////////////////////
2632/// Compute distance from array of input points having directions specified by dirs. Store output in dists
2633
2635 Double_t *step) const
2636{
2637 for (Int_t i = 0; i < vecsize; i++)
2638 dists[i] = DistFromInside(&points[3 * i], &dirs[3 * i], 3, step[i]);
2639}
2640
2641////////////////////////////////////////////////////////////////////////////////
2642/// Compute distance from array of input points having directions specified by dirs. Store output in dists
2643
2645 Double_t *step) const
2646{
2647 for (Int_t i = 0; i < vecsize; i++)
2648 dists[i] = DistFromOutside(&points[3 * i], &dirs[3 * i], 3, step[i]);
2649}
2650
2651////////////////////////////////////////////////////////////////////////////////
2652/// Compute safe distance from each of the points in the input array.
2653/// Input: Array of point coordinates, array of statuses for these points, size of the arrays
2654/// Output: Safety values
2655
2657{
2658 for (Int_t i = 0; i < vecsize; i++)
2659 safe[i] = Safety(&points[3 * i], inside[i]);
2660}
#define b(i)
Definition RSha256.hxx:100
#define c(i)
Definition RSha256.hxx:101
#define a(i)
Definition RSha256.hxx:99
#define s1(x)
Definition RSha256.hxx:91
int Int_t
Signed integer 4 bytes (int)
Definition RtypesCore.h:59
float Float_t
Float 4 bytes (float)
Definition RtypesCore.h:71
constexpr Bool_t kFALSE
Definition RtypesCore.h:108
double Double_t
Double 8 bytes.
Definition RtypesCore.h:73
constexpr Bool_t kTRUE
Definition RtypesCore.h:107
const char Option_t
Option string (const char)
Definition RtypesCore.h:80
ROOT::Detail::TRangeCast< T, true > TRangeDynCast
TRangeDynCast is an adapter class that allows the typed iteration through a TCollection.
void Error(const char *location, const char *msgfmt,...)
Use this function in case an error occurred.
Definition TError.cxx:208
void Fatal(const char *location, const char *msgfmt,...)
Use this function in case of a fatal error. It will abort the program.
Definition TError.cxx:267
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t r
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize id
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t points
char name[80]
Definition TGX11.cxx:110
R__EXTERN TGeoManager * gGeoManager
float xmin
float ymin
float xmax
float ymax
Generic 3D primitive description class.
Definition TBuffer3D.h:18
Int_t * fPols
Definition TBuffer3D.h:115
UInt_t NbPols() const
Definition TBuffer3D.h:82
UInt_t NbPnts() const
Definition TBuffer3D.h:80
UInt_t NbSegs() const
Definition TBuffer3D.h:81
Bool_t SectionsValid(UInt_t mask) const
Definition TBuffer3D.h:67
void SetSectionsValid(UInt_t mask)
Definition TBuffer3D.h:65
Int_t * fSegs
Definition TBuffer3D.h:114
Bool_t fLocalFrame
Definition TBuffer3D.h:90
Bool_t SetRawSizes(UInt_t reqPnts, UInt_t reqPntsCapacity, UInt_t reqSegs, UInt_t reqSegsCapacity, UInt_t reqPols, UInt_t reqPolsCapacity)
Set kRaw tessellation section of buffer with supplied sizes.
Double_t * fPnts
Definition TBuffer3D.h:113
Box class.
Definition TGeoBBox.h:17
void FillBuffer3D(TBuffer3D &buffer, Int_t reqSections, Bool_t localFrame) const override
Fills the supplied buffer, with sections in desired frame See TBuffer3D.h for explanation of sections...
Double_t fDX
Definition TGeoBBox.h:20
Double_t DistFromOutside(const Double_t *point, const Double_t *dir, Int_t iact=1, Double_t step=TGeoShape::Big(), Double_t *safe=nullptr) const override
Compute distance from outside point to surface of the box.
Definition TGeoBBox.cxx:432
Double_t fOrigin[3]
Definition TGeoBBox.h:23
void InspectShape() const override
Prints shape parameters.
Definition TGeoBBox.cxx:810
Double_t fDY
Definition TGeoBBox.h:21
Double_t fDZ
Definition TGeoBBox.h:22
A cone segment is a cone having a range in phi.
Definition TGeoCone.h:99
Double_t DistFromInside(const Double_t *point, const Double_t *dir, Int_t iact=1, Double_t step=TGeoShape::Big(), Double_t *safe=nullptr) const override
compute distance from inside point to surface of the tube segment
void InspectShape() const override
print shape parameters
TGeoShape * GetMakeRuntimeShape(TGeoShape *mother, TGeoMatrix *mat) const override
in case shape has some negative parameters, these has to be computed in order to fit the mother
static void ComputeNormalS(const Double_t *point, const Double_t *dir, Double_t *norm, Double_t dz, Double_t rmin1, Double_t rmax1, Double_t rmin2, Double_t rmax2, Double_t c1, Double_t s1, Double_t c2, Double_t s2)
Compute normal to closest surface from POINT.
Double_t fCdfi
Definition TGeoCone.h:111
Double_t Capacity() const override
Computes capacity of the shape in [length^3].
void Sizeof3D() const override
Fill size of this 3-D object.
Double_t fCm
Definition TGeoCone.h:110
static Double_t DistFromOutsideS(const Double_t *point, const Double_t *dir, Double_t dz, Double_t rmin1, Double_t rmax1, Double_t rmin2, Double_t rmax2, Double_t c1, Double_t s1, Double_t c2, Double_t s2, Double_t cm, Double_t sm, Double_t cdfi)
compute distance from outside point to surface of arbitrary tube
void SavePrimitive(std::ostream &out, Option_t *option="") override
Save a primitive as a C++ statement(s) on output stream "out".
TGeoConeSeg()
Default constructor.
Double_t GetAxisRange(Int_t iaxis, Double_t &xlo, Double_t &xhi) const override
Get range of shape for a given axis.
TGeoVolume * Divide(TGeoVolume *voldiv, const char *divname, Int_t iaxis, Int_t ndiv, Double_t start, Double_t step) override
Divide this cone segment shape belonging to volume "voldiv" into ndiv volumes called divname,...
void ComputeNormal(const Double_t *point, const Double_t *dir, Double_t *norm) const override
Compute normal to closest surface from POINT.
void ComputeBBox() override
compute bounding box of the tube segment
Bool_t Contains(const Double_t *point) const override
test if point is inside this sphere
Double_t fC1
Definition TGeoCone.h:106
void DistFromInside_v(const Double_t *points, const Double_t *dirs, Double_t *dists, Int_t vecsize, Double_t *step) const override
Compute distance from array of input points having directions specified by dirs. Store output in dist...
void SetSegsAndPols(TBuffer3D &buffer) const override
Fill TBuffer3D structure for segments and polygons.
void InitTrigonometry()
Init frequently used trigonometric values.
Int_t DistancetoPrimitive(Int_t px, Int_t py) override
compute closest distance from point px,py to each corner
void Contains_v(const Double_t *points, Bool_t *inside, Int_t vecsize) const override
Check the inside status for each of the points in the array.
void GetMeshNumbers(Int_t &nvert, Int_t &nsegs, Int_t &npols) const override
Returns numbers of vertices, segments and polygons composing the shape mesh.
static Double_t SafetyS(const Double_t *point, Bool_t in, Double_t dz, Double_t rmin1, Double_t rmax1, Double_t rmin2, Double_t rmax2, Double_t phi1, Double_t phi2, Int_t skipz=0)
Static method to compute the closest distance from given point to this shape.
Int_t GetNmeshVertices() const override
Return number of vertices of the mesh representation.
const TBuffer3D & GetBuffer3D(Int_t reqSections, Bool_t localFrame) const override
Fills a static 3D buffer and returns a reference.
void ComputeNormal_v(const Double_t *points, const Double_t *dirs, Double_t *norms, Int_t vecsize) override
Compute the normal for an array o points so that norm.dot.dir is positive Input: Arrays of point coor...
Double_t fPhi1
Definition TGeoCone.h:102
void SetConsDimensions(Double_t dz, Double_t rmin1, Double_t rmax1, Double_t rmin2, Double_t rmax2, Double_t phi1, Double_t phi2)
Set dimensions of the cone segment.
Double_t fS2
Definition TGeoCone.h:107
void AfterStreamer() override
Function called after streaming an object of this class.
void DistFromOutside_v(const Double_t *points, const Double_t *dirs, Double_t *dists, Int_t vecsize, Double_t *step) const override
Compute distance from array of input points having directions specified by dirs. Store output in dist...
~TGeoConeSeg() override
destructor
Double_t DistFromOutside(const Double_t *point, const Double_t *dir, Int_t iact=1, Double_t step=TGeoShape::Big(), Double_t *safe=nullptr) const override
compute distance from outside point to surface of the tube
Double_t fS1
Definition TGeoCone.h:105
static Double_t DistFromInsideS(const Double_t *point, const Double_t *dir, Double_t dz, Double_t rmin1, Double_t rmax1, Double_t rmin2, Double_t rmax2, Double_t c1, Double_t s1, Double_t c2, Double_t s2, Double_t cm, Double_t sm, Double_t cdfi)
compute distance from inside point to surface of the tube segment
static Double_t DistToCons(const Double_t *point, const Double_t *dir, Double_t r1, Double_t z1, Double_t r2, Double_t z2, Double_t phi1, Double_t phi2)
Static method to compute distance to a conical surface with :
Double_t Safety(const Double_t *point, Bool_t in=kTRUE) const override
computes the closest distance from given point to this shape, according to option.
Double_t fSm
Definition TGeoCone.h:109
Double_t fPhi2
Definition TGeoCone.h:103
void Safety_v(const Double_t *points, const Bool_t *inside, Double_t *safe, Int_t vecsize) const override
Compute safe distance from each of the points in the input array.
Double_t fC2
Definition TGeoCone.h:108
Bool_t GetPointsOnSegments(Int_t npoints, Double_t *array) const override
Fills array with n random points located on the line segments of the shape mesh.
void SetPoints(Double_t *points) const override
Create cone segment mesh points.
void SetDimensions(Double_t *param) override
Set dimensions of the cone segment from an array.
void GetBoundingCylinder(Double_t *param) const override
Fill vector param[4] with the bounding cylinder parameters.
TBuffer3D * MakeBuffer3D() const override
Creates a TBuffer3D describing this shape.
The cones are defined by 5 parameters:
Definition TGeoCone.h:17
void SetDimensions(Double_t *param) override
Set cone dimensions from an array.
Double_t Capacity() const override
Computes capacity of the shape in [length^3].
Definition TGeoCone.cxx:157
const char * GetAxisName(Int_t iaxis) const override
Returns name of axis IAXIS.
Definition TGeoCone.cxx:672
Double_t GetAxisRange(Int_t iaxis, Double_t &xlo, Double_t &xhi) const override
Get range of shape for a given axis.
Definition TGeoCone.cxx:685
static void ComputeNormalS(const Double_t *point, const Double_t *dir, Double_t *norm, Double_t dz, Double_t rmin1, Double_t rmax1, Double_t rmin2, Double_t rmax2)
Compute normal to closest surface from POINT.
Definition TGeoCone.cxx:233
void SetConeDimensions(Double_t dz, Double_t rmin1, Double_t rmax1, Double_t rmin2, Double_t rmax2)
Set cone dimensions.
Definition TGeoCone.cxx:993
static Double_t DistFromInsideS(const Double_t *point, const Double_t *dir, Double_t dz, Double_t rmin1, Double_t rmax1, Double_t rmin2, Double_t rmax2)
Compute distance from inside point to surface of the cone (static) Boundary safe algorithm.
Definition TGeoCone.cxx:286
~TGeoCone() override
destructor
Definition TGeoCone.cxx:175
void DistFromOutside_v(const Double_t *points, const Double_t *dirs, Double_t *dists, Int_t vecsize, Double_t *step) const override
Compute distance from array of input points having directions specified by dirs. Store output in dist...
Double_t fRmax1
Definition TGeoCone.h:22
Double_t fRmax2
Definition TGeoCone.h:24
void SetPoints(Double_t *points) const override
Create cone mesh points.
Int_t DistancetoPrimitive(Int_t px, Int_t py) override
compute closest distance from point px,py to each corner
Definition TGeoCone.cxx:602
Double_t Safety(const Double_t *point, Bool_t in=kTRUE) const override
computes the closest distance from given point to this shape, according to option.
Definition TGeoCone.cxx:922
Double_t DistFromOutside(const Double_t *point, const Double_t *dir, Int_t iact=1, Double_t step=TGeoShape::Big(), Double_t *safe=nullptr) const override
compute distance from outside point to surface of the tube
Definition TGeoCone.cxx:542
void DistFromInside_v(const Double_t *points, const Double_t *dirs, Double_t *dists, Int_t vecsize, Double_t *step) const override
Compute distance from array of input points having directions specified by dirs. Store output in dist...
void ComputeBBox() override
compute bounding box of the sphere
Definition TGeoCone.cxx:180
TGeoVolume * Divide(TGeoVolume *voldiv, const char *divname, Int_t iaxis, Int_t ndiv, Double_t start, Double_t step) override
Divide this cone shape belonging to volume "voldiv" into ndiv volumes called divname,...
Definition TGeoCone.cxx:618
static Double_t DistFromOutsideS(const Double_t *point, const Double_t *dir, Double_t dz, Double_t rmin1, Double_t rmax1, Double_t rmin2, Double_t rmax2)
Compute distance from outside point to surface of the tube Boundary safe algorithm.
Definition TGeoCone.cxx:388
Double_t fDz
Definition TGeoCone.h:20
void Contains_v(const Double_t *points, Bool_t *inside, Int_t vecsize) const override
Check the inside status for each of the points in the array.
TGeoCone()
Default constructor.
Definition TGeoCone.cxx:99
TGeoShape * GetMakeRuntimeShape(TGeoShape *mother, TGeoMatrix *mat) const override
in case shape has some negative parameters, these has to be computed in order to fit the mother
Definition TGeoCone.cxx:722
void GetMeshNumbers(Int_t &nvert, Int_t &nsegs, Int_t &npols) const override
Returns numbers of vertices, segments and polygons composing the shape mesh.
Bool_t GetPointsOnSegments(Int_t npoints, Double_t *array) const override
Fills array with n random points located on the line segments of the shape mesh.
Definition TGeoCone.cxx:755
Bool_t Contains(const Double_t *point) const override
test if point is inside this cone
Definition TGeoCone.cxx:270
static void DistToCone(const Double_t *point, const Double_t *dir, Double_t dz, Double_t r1, Double_t r2, Double_t &b, Double_t &delta)
Static method to compute distance to a conical surface with :
Definition TGeoCone.cxx:565
static Double_t SafetyS(const Double_t *point, Bool_t in, Double_t dz, Double_t rmin1, Double_t rmax1, Double_t rmin2, Double_t rmax2, Int_t skipz=0)
computes the closest distance from given point to this shape, according to option.
Definition TGeoCone.cxx:940
void Sizeof3D() const override
Fill size of this 3-D object.
void Safety_v(const Double_t *points, const Bool_t *inside, Double_t *safe, Int_t vecsize) const override
Compute safe distance from each of the points in the input array.
TBuffer3D * MakeBuffer3D() const override
Creates a TBuffer3D describing this shape.
Definition TGeoCone.cxx:823
void InspectShape() const override
print shape parameters
Definition TGeoCone.cxx:807
void ComputeNormal_v(const Double_t *points, const Double_t *dirs, Double_t *norms, Int_t vecsize) override
Compute the normal for an array o points so that norm.dot.dir is positive Input: Arrays of point coor...
void ComputeNormal(const Double_t *point, const Double_t *dir, Double_t *norm) const override
Compute normal to closest surface from POINT.
Definition TGeoCone.cxx:190
const TBuffer3D & GetBuffer3D(Int_t reqSections, Bool_t localFrame) const override
Fills a static 3D buffer and returns a reference.
Int_t GetNmeshVertices() const override
Return number of vertices of the mesh representation.
Double_t DistFromInside(const Double_t *point, const Double_t *dir, Int_t iact=1, Double_t step=TGeoShape::Big(), Double_t *safe=nullptr) const override
Compute distance from inside point to surface of the cone Boundary safe algorithm.
Definition TGeoCone.cxx:371
Double_t fRmin1
Definition TGeoCone.h:21
Double_t fRmin2
Definition TGeoCone.h:23
void SavePrimitive(std::ostream &out, Option_t *option="") override
Save a primitive as a C++ statement(s) on output stream "out".
Definition TGeoCone.cxx:975
void GetBoundingCylinder(Double_t *param) const override
Fill vector param[4] with the bounding cylinder parameters.
Definition TGeoCone.cxx:708
void SetSegsAndPols(TBuffer3D &buffer) const override
Fill TBuffer3D structure for segments and polygons.
Definition TGeoCone.cxx:842
TGeoVolumeMulti * MakeVolumeMulti(const char *name, TGeoMedium *medium)
Make a TGeoVolumeMulti handling a list of volumes.
Int_t GetNsegments() const
Get number of segments approximating circles.
Geometrical transformation package.
Definition TGeoMatrix.h:38
Node containing an offset.
Definition TGeoNode.h:184
a cylindrical phi divison pattern
base finder class for patterns. A pattern is specifying a division type
a Z axis divison pattern
Base abstract class for all shapes.
Definition TGeoShape.h:25
static Double_t Big()
Definition TGeoShape.h:94
Int_t GetBasicColor() const
Get the basic color (0-7).
void TransformPoints(Double_t *points, UInt_t NbPoints) const
Tranform a set of points (LocalToMaster)
void SetShapeBit(UInt_t f, Bool_t set)
Equivalent of TObject::SetBit.
static Double_t DistToPhiMin(const Double_t *point, const Double_t *dir, Double_t s1, Double_t c1, Double_t s2, Double_t c2, Double_t sm, Double_t cm, Bool_t in=kTRUE)
compute distance from point (inside phi) to both phi planes. Return minimum.
static Double_t SafetyPhi(const Double_t *point, Bool_t in, Double_t phi1, Double_t phi2)
Static method to compute safety w.r.t a phi corner defined by cosines/sines of the angles phi1,...
static Bool_t IsSameWithinTolerance(Double_t a, Double_t b)
Check if two numbers differ with less than a tolerance.
const char * GetPointerName() const
Provide a pointer name containing uid.
Int_t ShapeDistancetoPrimitive(Int_t numpoints, Int_t px, Int_t py) const
Returns distance to shape primitive mesh.
static void NormalPhi(const Double_t *point, const Double_t *dir, Double_t *norm, Double_t c1, Double_t s1, Double_t c2, Double_t s2)
Static method to compute normal to phi planes.
static Double_t SafetySeg(Double_t r, Double_t z, Double_t r1, Double_t z1, Double_t r2, Double_t z2, Bool_t outer)
Compute distance from point of coordinates (r,z) to segment (r1,z1):(r2,z2)
const char * GetName() const override
Get the shape name.
@ kGeoSavePrimitive
Definition TGeoShape.h:64
@ kGeoRunTimeShape
Definition TGeoShape.h:40
@ kGeoConeSeg
Definition TGeoShape.h:49
static Double_t Tolerance()
Definition TGeoShape.h:97
static Bool_t IsCloseToPhi(Double_t epsil, const Double_t *point, Double_t c1, Double_t s1, Double_t c2, Double_t s2)
True if point is closer than epsil to one of the phi planes defined by c1,s1 or c2,...
Bool_t TestShapeBit(UInt_t f) const
Definition TGeoShape.h:175
Volume families.
Definition TGeoVolume.h:266
TGeoVolume, TGeoVolumeMulti, TGeoVolumeAssembly are the volume classes.
Definition TGeoVolume.h:43
R__ALWAYS_INLINE Bool_t TestBit(UInt_t f) const
Definition TObject.h:202
virtual const char * ClassName() const
Returns name of class to which the object belongs.
Definition TObject.cxx:226
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition TObject.cxx:1057
void SetBit(UInt_t f, Bool_t set)
Set or unset the user status bits as specified in f.
Definition TObject.cxx:864
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition TObject.cxx:1071
Basic string class.
Definition TString.h:138
const char * Data() const
Definition TString.h:384
void box(Int_t pat, Double_t x1, Double_t y1, Double_t x2, Double_t y2)
Definition fillpatterns.C:1
return c1
Definition legend1.C:41
const Int_t n
Definition legend1.C:16
return c2
Definition legend2.C:14
Long64_t LocMin(Long64_t n, const T *a)
Returns index of array with the minimum element.
Definition TMath.h:993
Short_t Max(Short_t a, Short_t b)
Returns the largest of a and b.
Definition TMathBase.h:251
T1 Sign(T1 a, T2 b)
Returns a value with the magnitude of a and the sign of b.
Definition TMathBase.h:176
Double_t ATan2(Double_t y, Double_t x)
Returns the principal value of the arc tangent of y/x, expressed in radians.
Definition TMath.h:657
Long64_t LocMax(Long64_t n, const T *a)
Returns index of array with the maximum element.
Definition TMath.h:1099
constexpr Double_t DegToRad()
Conversion from degree to radian: .
Definition TMath.h:82
Double_t Sqrt(Double_t x)
Returns the square root of x.
Definition TMath.h:673
Short_t Min(Short_t a, Short_t b)
Returns the smallest of a and b.
Definition TMathBase.h:199
Double_t Cos(Double_t)
Returns the cosine of an angle of x radians.
Definition TMath.h:605
constexpr Double_t Pi()
Definition TMath.h:40
Double_t Sin(Double_t)
Returns the sine of an angle of x radians.
Definition TMath.h:599
constexpr Double_t RadToDeg()
Conversion from radian to degree: .
Definition TMath.h:75
Short_t Abs(Short_t d)
Returns the absolute value of parameter Short_t d.
Definition TMathBase.h:124
constexpr Double_t TwoPi()
Definition TMath.h:47