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