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