Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
TGeoPara.cxx
Go to the documentation of this file.
1 // @(#)root/geom:$Id$
2 // Author: Andrei Gheata 31/01/02
3 // TGeoPara::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 TGeoPara
14 \ingroup Geometry_classes
15 
16 Parallelepiped class. It has 6 parameters :
17 
18  - dx, dy, dz - half lengths in X, Y, Z
19  - alpha - angle w.r.t the Y axis from center of low Y edge to
20  center of high Y edge [deg]
21  - theta, phi - polar and azimuthal angles of the segment between
22  low and high Z surfaces [deg]
23 
24 Begin_Macro(source)
25 {
26  TCanvas *c = new TCanvas("c", "c",0,0,600,600);
27  new TGeoManager("para", "poza1");
28  TGeoMaterial *mat = new TGeoMaterial("Al", 26.98,13,2.7);
29  TGeoMedium *med = new TGeoMedium("MED",1,mat);
30  TGeoVolume *top = gGeoManager->MakeBox("TOP",med,100,100,100);
31  gGeoManager->SetTopVolume(top);
32  TGeoVolume *vol = gGeoManager->MakePara("PARA",med, 20,30,40,30,15,30);
33  vol->SetLineWidth(2);
34  top->AddNode(vol,1);
35  gGeoManager->CloseGeometry();
36  gGeoManager->SetNsegments(80);
37  top->Draw();
38  TView *view = gPad->GetView();
39  view->ShowAxis();
40 }
41 End_Macro
42 */
43 
44 #include "Riostream.h"
45 
46 #include "TGeoManager.h"
47 #include "TGeoMatrix.h"
48 #include "TGeoVolume.h"
49 #include "TGeoPara.h"
50 #include "TMath.h"
51 
52 ClassImp(TGeoPara);
53 
54 ////////////////////////////////////////////////////////////////////////////////
55 /// Default constructor
56 
57 TGeoPara::TGeoPara()
58 {
59  SetShapeBit(TGeoShape::kGeoPara);
60  fX = fY = fZ = 0;
61  fAlpha = 0;
62  fTheta = 0;
63  fPhi = 0;
64  fTxy = 0;
65  fTxz = 0;
66  fTyz = 0;
67 }
68 
69 ////////////////////////////////////////////////////////////////////////////////
70 /// Default constructor specifying minimum and maximum radius
71 
72 TGeoPara::TGeoPara(Double_t dx, Double_t dy, Double_t dz, Double_t alpha,
73  Double_t theta, Double_t phi)
74  :TGeoBBox(0, 0, 0)
75 {
76  SetShapeBit(TGeoShape::kGeoPara);
77  fX = dx;
78  fY = dy;
79  fZ = dz;
80  fAlpha = alpha;
81  fTheta = theta;
82  fPhi = phi;
83  fTxy = TMath::Tan(alpha*TMath::DegToRad());
84  Double_t tth = TMath::Tan(theta*TMath::DegToRad());
85  Double_t ph = phi*TMath::DegToRad();
86  fTxz = tth*TMath::Cos(ph);
87  fTyz = tth*TMath::Sin(ph);
88  if ((fX<0) || (fY<0) || (fZ<0)) {
89 // printf("para : %f %f %f\n", fX, fY, fZ);
90  SetShapeBit(kGeoRunTimeShape);
91  }
92  else ComputeBBox();
93 }
94 
95 ////////////////////////////////////////////////////////////////////////////////
96 /// Default constructor specifying minimum and maximum radius
97 
98 TGeoPara::TGeoPara(const char *name, Double_t dx, Double_t dy, Double_t dz, Double_t alpha,
99  Double_t theta, Double_t phi)
100  :TGeoBBox(name, 0, 0, 0)
101 {
102  SetShapeBit(TGeoShape::kGeoPara);
103  fX = dx;
104  fY = dy;
105  fZ = dz;
106  fAlpha = alpha;
107  fTheta = theta;
108  fPhi = phi;
109  fTxy = TMath::Tan(alpha*TMath::DegToRad());
110  Double_t tth = TMath::Tan(theta*TMath::DegToRad());
111  Double_t ph = phi*TMath::DegToRad();
112  fTxz = tth*TMath::Cos(ph);
113  fTyz = tth*TMath::Sin(ph);
114  if ((fX<0) || (fY<0) || (fZ<0)) {
115 // printf("para : %f %f %f\n", fX, fY, fZ);
116  SetShapeBit(kGeoRunTimeShape);
117  }
118  else ComputeBBox();
119 }
120 
121 ////////////////////////////////////////////////////////////////////////////////
122 /// Default constructor
123 /// - param[0] = dx
124 /// - param[1] = dy
125 /// - param[2] = dz
126 /// - param[3] = alpha
127 /// - param[4] = theta
128 /// - param[5] = phi
129 
130 TGeoPara::TGeoPara(Double_t *param)
131  :TGeoBBox(0, 0, 0)
132 {
133  SetShapeBit(TGeoShape::kGeoPara);
134  SetDimensions(param);
135  if ((fX<0) || (fY<0) || (fZ<0)) SetShapeBit(kGeoRunTimeShape);
136  else ComputeBBox();
137 }
138 
139 ////////////////////////////////////////////////////////////////////////////////
140 /// destructor
141 
142 TGeoPara::~TGeoPara()
143 {
144 }
145 
146 ////////////////////////////////////////////////////////////////////////////////
147 /// Computes capacity of the shape in [length^3]
148 
149 Double_t TGeoPara::Capacity() const
150 {
151  Double_t capacity = 8.*fX*fY*fZ;
152  return capacity;
153 }
154 
155 ////////////////////////////////////////////////////////////////////////////////
156 /// compute bounding box
157 
158 void TGeoPara::ComputeBBox()
159 {
160  Double_t dx = fX+fY*TMath::Abs(fTxy)+fZ*TMath::Abs(fTxz);
161  Double_t dy = fY+fZ*TMath::Abs(fTyz);
162  Double_t dz = fZ;
163  TGeoBBox::SetBoxDimensions(dx, dy, dz);
164  memset(fOrigin, 0, 3*sizeof(Double_t));
165 }
166 
167 ////////////////////////////////////////////////////////////////////////////////
168 /// Compute normal to closest surface from POINT.
169 
170 void TGeoPara::ComputeNormal(const Double_t *point, const Double_t *dir, Double_t *norm)
171 {
172  Double_t saf[3];
173  // distance from point to higher Z face
174  saf[0] = TMath::Abs(fZ-TMath::Abs(point[2])); // Z
175 
176  Double_t yt = point[1]-fTyz*point[2];
177  saf[1] = TMath::Abs(fY-TMath::Abs(yt)); // Y
178  // cos of angle YZ
179  Double_t cty = 1.0/TMath::Sqrt(1.0+fTyz*fTyz);
180 
181  Double_t xt = point[0]-fTxz*point[2]-fTxy*yt;
182  saf[2] = TMath::Abs(fX-TMath::Abs(xt)); // X
183  // cos of angle XZ
184  Double_t ctx = 1.0/TMath::Sqrt(1.0+fTxy*fTxy+fTxz*fTxz);
185  saf[2] *= ctx;
186  saf[1] *= cty;
187  Int_t i = TMath::LocMin(3,saf);
188  switch (i) {
189  case 0:
190  norm[0] = norm[1] = 0;
191  norm[2] = TMath::Sign(1.,dir[2]);
192  return;
193  case 1:
194  norm[0] = 0;
195  norm[1] = cty;
196  norm[2] = - fTyz*cty;
197  break;
198  case 2:
199  norm[0] = TMath::Cos(fTheta*TMath::DegToRad())*TMath::Cos(fAlpha*TMath::DegToRad());
200  norm[1] = - TMath::Cos(fTheta*TMath::DegToRad())*TMath::Sin(fAlpha*TMath::DegToRad());
201  norm[2] = -TMath::Sin(fTheta*TMath::DegToRad());
202  }
203  if (norm[0]*dir[0]+norm[1]*dir[1]+norm[2]*dir[2]<0) {
204  norm[0] = -norm[0];
205  norm[1] = -norm[1];
206  norm[2] = -norm[2];
207  }
208 }
209 
210 ////////////////////////////////////////////////////////////////////////////////
211 /// test if point is inside this sphere
212 /// test Z range
213 
214 Bool_t TGeoPara::Contains(const Double_t *point) const
215 {
216  if (TMath::Abs(point[2]) > fZ) return kFALSE;
217  // check X and Y
218  Double_t yt=point[1]-fTyz*point[2];
219  if (TMath::Abs(yt) > fY) return kFALSE;
220  Double_t xt=point[0]-fTxz*point[2]-fTxy*yt;
221  if (TMath::Abs(xt) > fX) return kFALSE;
222  return kTRUE;
223 }
224 
225 ////////////////////////////////////////////////////////////////////////////////
226 /// compute distance from inside point to surface of the para
227 /// Boundary safe algorithm.
228 
229 Double_t TGeoPara::DistFromInside(const Double_t *point, const Double_t *dir, Int_t iact, Double_t step, Double_t *safe) const
230 {
231  if (iact<3 && safe) {
232  // compute safety
233  *safe = Safety(point, kTRUE);
234  if (iact==0) return TGeoShape::Big();
235  if (iact==1 && step<*safe) return TGeoShape::Big();
236  }
237  Double_t saf[2];
238  Double_t snxt = TGeoShape::Big();
239  Double_t s;
240  saf[0] = fZ+point[2];
241  saf[1] = fZ-point[2];
242  if (!TGeoShape::IsSameWithinTolerance(dir[2],0)) {
243  s = (dir[2]>0)?(saf[1]/dir[2]):(-saf[0]/dir[2]);
244  if (s<0) return 0.0;
245  if (s<snxt) snxt = s;
246  }
247  // distance from point to center axis on Y
248  Double_t yt = point[1]-fTyz*point[2];
249  saf[0] = fY+yt;
250  saf[1] = fY-yt;
251  Double_t dy = dir[1]-fTyz*dir[2];
252  if (!TGeoShape::IsSameWithinTolerance(dy,0)) {
253  s = (dy>0)?(saf[1]/dy):(-saf[0]/dy);
254  if (s<0) return 0.0;
255  if (s<snxt) snxt = s;
256  }
257  // distance from point to center axis on X
258  Double_t xt = point[0]-fTxz*point[2]-fTxy*yt;
259  saf[0] = fX+xt;
260  saf[1] = fX-xt;
261  Double_t dx = dir[0]-fTxz*dir[2]-fTxy*dy;
262  if (!TGeoShape::IsSameWithinTolerance(dx,0)) {
263  s = (dx>0)?(saf[1]/dx):(-saf[0]/dx);
264  if (s<0) return 0.0;
265  if (s<snxt) snxt = s;
266  }
267  return snxt;
268 }
269 
270 ////////////////////////////////////////////////////////////////////////////////
271 /// compute distance from inside point to surface of the para
272 
273 Double_t TGeoPara::DistFromOutside(const Double_t *point, const Double_t *dir, Int_t iact, Double_t step, Double_t *safe) const
274 {
275  Double_t snxt=TGeoShape::Big();
276  if (iact<3 && safe) {
277  // compute safe distance
278  *safe = Safety(point, kFALSE);
279  if (iact==0) return TGeoShape::Big();
280  if (iact==1 && step<*safe) return TGeoShape::Big();
281  }
282  Bool_t in = kTRUE;
283  Double_t safz;
284  safz = TMath::Abs(point[2])-fZ;
285  if (safz>0) {
286  // outside Z
287  if (point[2]*dir[2]>=0) return TGeoShape::Big();
288  in = kFALSE;
289  }
290  Double_t yt=point[1]-fTyz*point[2];
291  Double_t safy = TMath::Abs(yt)-fY;
292  Double_t dy=dir[1]-fTyz*dir[2];
293  if (safy>0) {
294  if (yt*dy>=0) return TGeoShape::Big();
295  in = kFALSE;
296  }
297  Double_t xt=point[0]-fTxy*yt-fTxz*point[2];
298  Double_t safx = TMath::Abs(xt)-fX;
299  Double_t dx=dir[0]-fTxy*dy-fTxz*dir[2];
300  if (safx>0) {
301  if (xt*dx>=0) return TGeoShape::Big();
302  in = kFALSE;
303  }
304  // protection in case point is actually inside
305  if (in) {
306  if (safz>safx && safz>safy) {
307  if (point[2]*dir[2]>0) return TGeoShape::Big();
308  return 0.0;
309  }
310  if (safx>safy) {
311  if (xt*dx>0) return TGeoShape::Big();
312  return 0.0;
313  }
314  if (yt*dy>0) return TGeoShape::Big();
315  return 0.0;
316  }
317  Double_t xnew,ynew,znew;
318  if (safz>0) {
319  snxt = safz/TMath::Abs(dir[2]);
320  xnew = point[0]+snxt*dir[0];
321  ynew = point[1]+snxt*dir[1];
322  znew = (point[2]>0)?fZ:(-fZ);
323  Double_t ytn = ynew-fTyz*znew;
324  if (TMath::Abs(ytn)<=fY) {
325  Double_t xtn = xnew-fTxy*ytn-fTxz*znew;
326  if (TMath::Abs(xtn)<=fX) return snxt;
327  }
328  }
329  if (safy>0) {
330  snxt = safy/TMath::Abs(dy);
331  znew = point[2]+snxt*dir[2];
332  if (TMath::Abs(znew)<=fZ) {
333  Double_t ytn = (yt>0)?fY:(-fY);
334  xnew = point[0]+snxt*dir[0];
335  Double_t xtn = xnew-fTxy*ytn-fTxz*znew;
336  if (TMath::Abs(xtn)<=fX) return snxt;
337  }
338  }
339  if (safx>0) {
340  snxt = safx/TMath::Abs(dx);
341  znew = point[2]+snxt*dir[2];
342  if (TMath::Abs(znew)<=fZ) {
343  ynew = point[1]+snxt*dir[1];
344  Double_t ytn = ynew-fTyz*znew;
345  if (TMath::Abs(ytn)<=fY) return snxt;
346  }
347  }
348  return TGeoShape::Big();
349 }
350 
351 ////////////////////////////////////////////////////////////////////////////////
352 /// Divide this parallelepiped shape belonging to volume "voldiv" into ndiv equal volumes
353 /// called divname, from start position with the given step. Returns pointer
354 /// to created division cell volume. In case a wrong division axis is supplied,
355 /// returns pointer to volume to be divided.
356 
357 TGeoVolume *TGeoPara::Divide(TGeoVolume *voldiv, const char *divname, Int_t iaxis, Int_t ndiv,
358  Double_t start, Double_t step)
359 {
360  TGeoShape *shape; //--- shape to be created
361  TGeoVolume *vol; //--- division volume to be created
362  TGeoVolumeMulti *vmulti; //--- generic divided volume
363  TGeoPatternFinder *finder; //--- finder to be attached
364  TString opt = ""; //--- option to be attached
365  Double_t end=start+ndiv*step;
366  switch (iaxis) {
367  case 1: //--- divide on X
368  shape = new TGeoPara(step/2, fY, fZ,fAlpha,fTheta, fPhi);
369  finder = new TGeoPatternParaX(voldiv, ndiv, start, end);
370  opt = "X";
371  break;
372  case 2: //--- divide on Y
373  shape = new TGeoPara(fX, step/2, fZ, fAlpha, fTheta, fPhi);
374  finder = new TGeoPatternParaY(voldiv, ndiv, start, end);
375  opt = "Y";
376  break;
377  case 3: //--- divide on Z
378  shape = new TGeoPara(fX, fY, step/2, fAlpha, fTheta, fPhi);
379  finder = new TGeoPatternParaZ(voldiv, ndiv, start, end);
380  opt = "Z";
381  break;
382  default:
383  Error("Divide", "Wrong axis type for division");
384  return 0;
385  }
386  vol = new TGeoVolume(divname, shape, voldiv->GetMedium());
387  vmulti = gGeoManager->MakeVolumeMulti(divname, voldiv->GetMedium());
388  vmulti->AddVolume(vol);
389  voldiv->SetFinder(finder);
390  finder->SetDivIndex(voldiv->GetNdaughters());
391  for (Int_t ic=0; ic<ndiv; ic++) {
392  voldiv->AddNodeOffset(vol, ic, start+step/2.+ic*step, opt.Data());
393  ((TGeoNodeOffset*)voldiv->GetNodes()->At(voldiv->GetNdaughters()-1))->SetFinder(finder);
394  }
395  return vmulti;
396 }
397 
398 ////////////////////////////////////////////////////////////////////////////////
399 /// Get range of shape for a given axis.
400 
401 Double_t TGeoPara::GetAxisRange(Int_t iaxis, Double_t &xlo, Double_t &xhi) const
402 {
403  xlo = 0;
404  xhi = 0;
405  Double_t dx = 0;
406  switch (iaxis) {
407  case 1:
408  xlo = -fX;
409  xhi = fX;
410  dx = xhi-xlo;
411  return dx;
412  case 2:
413  xlo = -fY;
414  xhi = fY;
415  dx = xhi-xlo;
416  return dx;
417  case 3:
418  xlo = -fZ;
419  xhi = fZ;
420  dx = xhi-xlo;
421  return dx;
422  }
423  return dx;
424 }
425 
426 ////////////////////////////////////////////////////////////////////////////////
427 /// Fill vector param[4] with the bounding cylinder parameters. The order
428 /// is the following : Rmin, Rmax, Phi1, Phi2
429 
430 void TGeoPara::GetBoundingCylinder(Double_t *param) const
431 {
432  TGeoBBox::GetBoundingCylinder(param);
433 }
434 
435 ////////////////////////////////////////////////////////////////////////////////
436 /// Fills real parameters of a positioned box inside this. Returns 0 if successful.
437 
438 Int_t TGeoPara::GetFittingBox(const TGeoBBox *parambox, TGeoMatrix *mat, Double_t &dx, Double_t &dy, Double_t &dz) const
439 {
440  dx=dy=dz=0;
441  if (mat->IsRotation()) {
442  Error("GetFittingBox", "cannot handle parametrized rotated volumes");
443  return 1; // ### rotation not accepted ###
444  }
445  //--> translate the origin of the parametrized box to the frame of this box.
446  Double_t origin[3];
447  mat->LocalToMaster(parambox->GetOrigin(), origin);
448  if (!Contains(origin)) {
449  Error("GetFittingBox", "wrong matrix - parametrized box is outside this");
450  return 1; // ### wrong matrix ###
451  }
452  //--> now we have to get the valid range for all parametrized axis
453  Double_t dd[3];
454  dd[0] = parambox->GetDX();
455  dd[1] = parambox->GetDY();
456  dd[2] = parambox->GetDZ();
457  //-> check if Z range is fixed
458  if (dd[2]<0) {
459  dd[2] = TMath::Min(origin[2]+fZ, fZ-origin[2]);
460  if (dd[2]<0) {
461  Error("GetFittingBox", "wrong matrix");
462  return 1;
463  }
464  }
465  if (dd[0]>=0 && dd[1]>=0) {
466  dx = dd[0];
467  dy = dd[1];
468  dz = dd[2];
469  return 0;
470  }
471  //-> check now range at Z = origin[2] +/- dd[2]
472  Double_t upper[8];
473  Double_t lower[8];
474  Double_t z=origin[2]-dd[2];
475  lower[0]=z*fTxz-fTxy*fY-fX;
476  lower[1]=-fY+z*fTyz;
477  lower[2]=z*fTxz+fTxy*fY-fX;
478  lower[3]=fY+z*fTyz;
479  lower[4]=z*fTxz+fTxy*fY+fX;
480  lower[5]=fY+z*fTyz;
481  lower[6]=z*fTxz-fTxy*fY+fX;
482  lower[7]=-fY+z*fTyz;
483  z=origin[2]+dd[2];
484  upper[0]=z*fTxz-fTxy*fY-fX;
485  upper[1]=-fY+z*fTyz;
486  upper[2]=z*fTxz+fTxy*fY-fX;
487  upper[3]=fY+z*fTyz;
488  upper[4]=z*fTxz+fTxy*fY+fX;
489  upper[5]=fY+z*fTyz;
490  upper[6]=z*fTxz-fTxy*fY+fX;
491  upper[7]=-fY+z*fTyz;
492 
493  Double_t ddmin=TGeoShape::Big();
494  for (Int_t iaxis=0; iaxis<2; iaxis++) {
495  if (dd[iaxis]>=0) continue;
496  ddmin=TGeoShape::Big();
497  for (Int_t ivert=0; ivert<4; ivert++) {
498  ddmin = TMath::Min(ddmin, TMath::Abs(origin[iaxis]-lower[2*ivert+iaxis]));
499  ddmin = TMath::Min(ddmin, TMath::Abs(origin[iaxis]-upper[2*ivert+iaxis]));
500  }
501  dd[iaxis] = ddmin;
502  }
503  dx = dd[0];
504  dy = dd[1];
505  dz = dd[2];
506  return 0;
507 }
508 
509 ////////////////////////////////////////////////////////////////////////////////
510 /// in case shape has some negative parameters, these has to be computed
511 /// in order to fit the mother
512 
513 TGeoShape *TGeoPara::GetMakeRuntimeShape(TGeoShape *mother, TGeoMatrix * /*mat*/) const
514 {
515  if (!TestShapeBit(kGeoRunTimeShape)) return 0;
516  if (!mother->TestShapeBit(kGeoPara)) {
517  Error("GetMakeRuntimeShape", "invalid mother");
518  return 0;
519  }
520  Double_t dx, dy, dz;
521  if (fX<0) dx=((TGeoPara*)mother)->GetX();
522  else dx=fX;
523  if (fY<0) dy=((TGeoPara*)mother)->GetY();
524  else dy=fY;
525  if (fZ<0) dz=((TGeoPara*)mother)->GetZ();
526  else dz=fZ;
527  return (new TGeoPara(dx, dy, dz, fAlpha, fTheta, fPhi));
528 }
529 
530 ////////////////////////////////////////////////////////////////////////////////
531 /// print shape parameters
532 
533 void TGeoPara::InspectShape() const
534 {
535  printf("*** Shape %s: TGeoPara ***\n", GetName());
536  printf(" dX = %11.5f\n", fX);
537  printf(" dY = %11.5f\n", fY);
538  printf(" dZ = %11.5f\n", fZ);
539  printf(" alpha = %11.5f\n", fAlpha);
540  printf(" theta = %11.5f\n", fTheta);
541  printf(" phi = %11.5f\n", fPhi);
542  printf(" Bounding box:\n");
543  TGeoBBox::InspectShape();
544 }
545 
546 ////////////////////////////////////////////////////////////////////////////////
547 /// computes the closest distance from given point to this shape, according
548 /// to option. The matching point on the shape is stored in spoint.
549 
550 Double_t TGeoPara::Safety(const Double_t *point, Bool_t in) const
551 {
552  Double_t saf[3];
553  // distance from point to higher Z face
554  saf[0] = fZ-TMath::Abs(point[2]); // Z
555 
556  Double_t yt = point[1]-fTyz*point[2];
557  saf[1] = fY-TMath::Abs(yt); // Y
558  // cos of angle YZ
559  Double_t cty = 1.0/TMath::Sqrt(1.0+fTyz*fTyz);
560 
561  Double_t xt = point[0]-fTxz*point[2]-fTxy*yt;
562  saf[2] = fX-TMath::Abs(xt); // X
563  // cos of angle XZ
564  Double_t ctx = 1.0/TMath::Sqrt(1.0+fTxy*fTxy+fTxz*fTxz);
565  saf[2] *= ctx;
566  saf[1] *= cty;
567  if (in) return saf[TMath::LocMin(3,saf)];
568  for (Int_t i=0; i<3; i++) saf[i]=-saf[i];
569  return saf[TMath::LocMax(3,saf)];
570 }
571 
572 ////////////////////////////////////////////////////////////////////////////////
573 /// Save a primitive as a C++ statement(s) on output stream "out".
574 
575 void TGeoPara::SavePrimitive(std::ostream &out, Option_t * /*option*/ /*= ""*/)
576 {
577  if (TObject::TestBit(kGeoSavePrimitive)) return;
578  out << " // Shape: " << GetName() << " type: " << ClassName() << std::endl;
579  out << " dx = " << fX << ";" << std::endl;
580  out << " dy = " << fY << ";" << std::endl;
581  out << " dz = " << fZ << ";" << std::endl;
582  out << " alpha = " << fAlpha<< ";" << std::endl;
583  out << " theta = " << fTheta << ";" << std::endl;
584  out << " phi = " << fPhi << ";" << std::endl;
585  out << " TGeoShape *" << GetPointerName() << " = new TGeoPara(\"" << GetName() << "\",dx,dy,dz,alpha,theta,phi);" << std::endl;
586  TObject::SetBit(TGeoShape::kGeoSavePrimitive);
587 }
588 
589 ////////////////////////////////////////////////////////////////////////////////
590 /// Set dimensions starting from an array.
591 
592 void TGeoPara::SetDimensions(Double_t *param)
593 {
594  fX = param[0];
595  fY = param[1];
596  fZ = param[2];
597  fAlpha = param[3];
598  fTheta = param[4];
599  fPhi = param[5];
600  fTxy = TMath::Tan(param[3]*TMath::DegToRad());
601  Double_t tth = TMath::Tan(param[4]*TMath::DegToRad());
602  Double_t ph = param[5]*TMath::DegToRad();
603  fTxz = tth*TMath::Cos(ph);
604  fTyz = tth*TMath::Sin(ph);
605 }
606 
607 ////////////////////////////////////////////////////////////////////////////////
608 /// Create PARA mesh points
609 
610 void TGeoPara::SetPoints(Double_t *points) const
611 {
612  if (!points) return;
613  Double_t txy = fTxy;
614  Double_t txz = fTxz;
615  Double_t tyz = fTyz;
616  *points++ = -fZ*txz-txy*fY-fX; *points++ = -fY-fZ*tyz; *points++ = -fZ;
617  *points++ = -fZ*txz+txy*fY-fX; *points++ = +fY-fZ*tyz; *points++ = -fZ;
618  *points++ = -fZ*txz+txy*fY+fX; *points++ = +fY-fZ*tyz; *points++ = -fZ;
619  *points++ = -fZ*txz-txy*fY+fX; *points++ = -fY-fZ*tyz; *points++ = -fZ;
620  *points++ = +fZ*txz-txy*fY-fX; *points++ = -fY+fZ*tyz; *points++ = +fZ;
621  *points++ = +fZ*txz+txy*fY-fX; *points++ = +fY+fZ*tyz; *points++ = +fZ;
622  *points++ = +fZ*txz+txy*fY+fX; *points++ = +fY+fZ*tyz; *points++ = +fZ;
623  *points++ = +fZ*txz-txy*fY+fX; *points++ = -fY+fZ*tyz; *points++ = +fZ;
624 }
625 
626 ////////////////////////////////////////////////////////////////////////////////
627 /// create sphere mesh points
628 
629 void TGeoPara::SetPoints(Float_t *points) const
630 {
631  if (!points) return;
632  Double_t txy = fTxy;
633  Double_t txz = fTxz;
634  Double_t tyz = fTyz;
635  *points++ = -fZ*txz-txy*fY-fX; *points++ = -fY-fZ*tyz; *points++ = -fZ;
636  *points++ = -fZ*txz+txy*fY-fX; *points++ = +fY-fZ*tyz; *points++ = -fZ;
637  *points++ = -fZ*txz+txy*fY+fX; *points++ = +fY-fZ*tyz; *points++ = -fZ;
638  *points++ = -fZ*txz-txy*fY+fX; *points++ = -fY-fZ*tyz; *points++ = -fZ;
639  *points++ = +fZ*txz-txy*fY-fX; *points++ = -fY+fZ*tyz; *points++ = +fZ;
640  *points++ = +fZ*txz+txy*fY-fX; *points++ = +fY+fZ*tyz; *points++ = +fZ;
641  *points++ = +fZ*txz+txy*fY+fX; *points++ = +fY+fZ*tyz; *points++ = +fZ;
642  *points++ = +fZ*txz-txy*fY+fX; *points++ = -fY+fZ*tyz; *points++ = +fZ;
643 }
644 
645 ////////////////////////////////////////////////////////////////////////////////
646 /// fill size of this 3-D object
647 
648 void TGeoPara::Sizeof3D() const
649 {
650  TGeoBBox::Sizeof3D();
651 }
652 
653 ////////////////////////////////////////////////////////////////////////////////
654 /// Check the inside status for each of the points in the array.
655 /// Input: Array of point coordinates + vector size
656 /// Output: Array of Booleans for the inside of each point
657 
658 void TGeoPara::Contains_v(const Double_t *points, Bool_t *inside, Int_t vecsize) const
659 {
660  for (Int_t i=0; i<vecsize; i++) inside[i] = Contains(&points[3*i]);
661 }
662 
663 ////////////////////////////////////////////////////////////////////////////////
664 /// Compute the normal for an array o points so that norm.dot.dir is positive
665 /// Input: Arrays of point coordinates and directions + vector size
666 /// Output: Array of normal directions
667 
668 void TGeoPara::ComputeNormal_v(const Double_t *points, const Double_t *dirs, Double_t *norms, Int_t vecsize)
669 {
670  for (Int_t i=0; i<vecsize; i++) ComputeNormal(&points[3*i], &dirs[3*i], &norms[3*i]);
671 }
672 
673 ////////////////////////////////////////////////////////////////////////////////
674 /// Compute distance from array of input points having directions specified by dirs. Store output in dists
675 
676 void TGeoPara::DistFromInside_v(const Double_t *points, const Double_t *dirs, Double_t *dists, Int_t vecsize, Double_t* step) const
677 {
678  for (Int_t i=0; i<vecsize; i++) dists[i] = DistFromInside(&points[3*i], &dirs[3*i], 3, step[i]);
679 }
680 
681 ////////////////////////////////////////////////////////////////////////////////
682 /// Compute distance from array of input points having directions specified by dirs. Store output in dists
683 
684 void TGeoPara::DistFromOutside_v(const Double_t *points, const Double_t *dirs, Double_t *dists, Int_t vecsize, Double_t* step) const
685 {
686  for (Int_t i=0; i<vecsize; i++) dists[i] = DistFromOutside(&points[3*i], &dirs[3*i], 3, step[i]);
687 }
688 
689 ////////////////////////////////////////////////////////////////////////////////
690 /// Compute safe distance from each of the points in the input array.
691 /// Input: Array of point coordinates, array of statuses for these points, size of the arrays
692 /// Output: Safety values
693 
694 void TGeoPara::Safety_v(const Double_t *points, const Bool_t *inside, Double_t *safe, Int_t vecsize) const
695 {
696  for (Int_t i=0; i<vecsize; i++) safe[i] = Safety(&points[3*i], inside[i]);
697 }