Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
TGeoTrd2.cxx
Go to the documentation of this file.
1 // @(#)root/geom:$Id$
2 // Author: Andrei Gheata 31/01/02
3 // TGeoTrd2::Contains() and DistFromInside() implemented by Mihaela Gheata
4 
5 /*************************************************************************
6  * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers. *
7  * All rights reserved. *
8  * *
9  * For the licensing terms see $ROOTSYS/LICENSE. *
10  * For the list of contributors see $ROOTSYS/README/CREDITS. *
11  *************************************************************************/
12 
13 /** \class TGeoTrd2
14 \ingroup Geometry_classes
15 A trapezoid with both x and y lengths varying with z. It
16 has 5 parameters, the half lengths in x at -dz and +dz, the half
17 lengths in y at -dz and +dz, and the half length in z (dz).
18 
19 Begin_Macro(source)
20 {
21  TCanvas *c = new TCanvas("c", "c",0,0,600,600);
22  new TGeoManager("trd2", "poza9");
23  TGeoMaterial *mat = new TGeoMaterial("Al", 26.98,13,2.7);
24  TGeoMedium *med = new TGeoMedium("MED",1,mat);
25  TGeoVolume *top = gGeoManager->MakeBox("TOP",med,100,100,100);
26  gGeoManager->SetTopVolume(top);
27  TGeoVolume *vol = gGeoManager->MakeTrd2("Trd2",med, 10,20,30,10,40);
28  vol->SetLineWidth(2);
29  gGeoManager->CloseGeometry();
30  gGeoManager->SetNsegments(80);
31  top->Draw();
32  TView *view = gPad->GetView();
33  view->ShowAxis();
34 }
35 End_Macro
36 */
37 
38 #include "Riostream.h"
39 
40 #include "TGeoManager.h"
41 #include "TGeoMatrix.h"
42 #include "TGeoVolume.h"
43 #include "TGeoTrd2.h"
44 #include "TMath.h"
45 
46 ClassImp(TGeoTrd2);
47 
48 ////////////////////////////////////////////////////////////////////////////////
49 /// dummy ctor
50 
51 TGeoTrd2::TGeoTrd2()
52 {
53  SetShapeBit(kGeoTrd2);
54  fDz = fDx1 = fDx2 = fDy1 = fDy2 = 0;
55 }
56 
57 ////////////////////////////////////////////////////////////////////////////////
58 /// constructor.
59 
60 TGeoTrd2::TGeoTrd2(Double_t dx1, Double_t dx2, Double_t dy1, Double_t dy2, Double_t dz)
61  :TGeoBBox(0,0,0)
62 {
63  SetShapeBit(kGeoTrd2);
64  fDx1 = dx1;
65  fDx2 = dx2;
66  fDy1 = dy1;
67  fDy2 = dy2;
68  fDz = dz;
69  if ((fDx1<0) || (fDx2<0) || (fDy1<0) || (fDy2<0) || (fDz<0)) {
70  SetShapeBit(kGeoRunTimeShape);
71  printf("trd2 : dx1=%f, dx2=%f, dy1=%f, dy2=%f, dz=%f\n",
72  dx1,dx2,dy1,dy2,dz);
73  }
74  else ComputeBBox();
75 }
76 
77 ////////////////////////////////////////////////////////////////////////////////
78 /// constructor.
79 
80 TGeoTrd2::TGeoTrd2(const char * name, Double_t dx1, Double_t dx2, Double_t dy1, Double_t dy2, Double_t dz)
81  :TGeoBBox(name, 0,0,0)
82 {
83  SetShapeBit(kGeoTrd2);
84  fDx1 = dx1;
85  fDx2 = dx2;
86  fDy1 = dy1;
87  fDy2 = dy2;
88  fDz = dz;
89  if ((fDx1<0) || (fDx2<0) || (fDy1<0) || (fDy2<0) || (fDz<0)) {
90  SetShapeBit(kGeoRunTimeShape);
91  printf("trd2 : dx1=%f, dx2=%f, dy1=%f, dy2=%f, dz=%f\n",
92  dx1,dx2,dy1,dy2,dz);
93  }
94  else ComputeBBox();
95 }
96 
97 ////////////////////////////////////////////////////////////////////////////////
98 /// ctor with an array of parameters
99 /// - param[0] = dx1
100 /// - param[1] = dx2
101 /// - param[2] = dy1
102 /// - param[3] = dy2
103 /// - param[4] = dz
104 
105 TGeoTrd2::TGeoTrd2(Double_t *param)
106  :TGeoBBox(0,0,0)
107 {
108  SetShapeBit(kGeoTrd2);
109  SetDimensions(param);
110  if ((fDx1<0) || (fDx2<0) || (fDy1<0) || (fDy2<0) || (fDz<0)) SetShapeBit(kGeoRunTimeShape);
111  else ComputeBBox();
112 }
113 
114 ////////////////////////////////////////////////////////////////////////////////
115 /// destructor
116 
117 TGeoTrd2::~TGeoTrd2()
118 {
119 }
120 
121 ////////////////////////////////////////////////////////////////////////////////
122 /// Computes capacity of the shape in [length^3]
123 
124 Double_t TGeoTrd2::Capacity() const
125 {
126  Double_t capacity = 2*(fDx1+fDx2)*(fDy1+fDy2)*fDz +
127  (2./3.)*(fDx1-fDx2)*(fDy1-fDy2)*fDz;
128  return capacity;
129 }
130 
131 ////////////////////////////////////////////////////////////////////////////////
132 /// compute bounding box for a trd2
133 
134 void TGeoTrd2::ComputeBBox()
135 {
136  fDX = TMath::Max(fDx1, fDx2);
137  fDY = TMath::Max(fDy1, fDy2);
138  fDZ = fDz;
139  memset(fOrigin, 0, 3*sizeof(Double_t));
140 }
141 
142 ////////////////////////////////////////////////////////////////////////////////
143 /// Compute normal to closest surface from POINT.
144 
145 void TGeoTrd2::ComputeNormal(const Double_t *point, const Double_t *dir, Double_t *norm)
146 {
147  Double_t safe, safemin;
148  Double_t fx = 0.5*(fDx1-fDx2)/fDz;
149  Double_t calf = 1./TMath::Sqrt(1.0+fx*fx);
150  // check Z facettes
151  safe = safemin = TMath::Abs(fDz-TMath::Abs(point[2]));
152  norm[0] = norm[1] = 0;
153  norm[2] = (dir[2]>=0)?1:-1;
154  if (safe<TGeoShape::Tolerance()) return;
155  // check X facettes
156  Double_t distx = 0.5*(fDx1+fDx2)-fx*point[2];
157  if (distx>=0) {
158  safe=TMath::Abs(distx-TMath::Abs(point[0]))*calf;
159  if (safe<safemin) {
160  safemin = safe;
161  norm[0] = (point[0]>0)?calf:(-calf);
162  norm[1] = 0;
163  norm[2] = calf*fx;
164  Double_t dot = norm[0]*dir[0]+norm[1]*dir[1]+norm[2]*dir[2];
165  if (dot<0) {
166  norm[0]=-norm[0];
167  norm[2]=-norm[2];
168  }
169  if (safe<TGeoShape::Tolerance()) return;
170  }
171  }
172 
173  Double_t fy = 0.5*(fDy1-fDy2)/fDz;
174  calf = 1./TMath::Sqrt(1.0+fy*fy);
175 
176  // check Y facettes
177  distx = 0.5*(fDy1+fDy2)-fy*point[2];
178  if (distx>=0) {
179  safe=TMath::Abs(distx-TMath::Abs(point[1]))*calf;
180  if (safe<safemin) {
181  norm[0] = 0;
182  norm[1] = (point[1]>0)?calf:(-calf);
183  norm[2] = calf*fy;
184  Double_t dot = norm[0]*dir[0]+norm[1]*dir[1]+norm[2]*dir[2];
185  if (dot<0) {
186  norm[1]=-norm[1];
187  norm[2]=-norm[2];
188  }
189  }
190  }
191 }
192 
193 ////////////////////////////////////////////////////////////////////////////////
194 /// test if point is inside this shape
195 /// check Z range
196 
197 Bool_t TGeoTrd2::Contains(const Double_t *point) const
198 {
199  if (TMath::Abs(point[2]) > fDz) return kFALSE;
200  // then y
201  Double_t dy = 0.5*(fDy2*(point[2]+fDz)+fDy1*(fDz-point[2]))/fDz;
202  if (TMath::Abs(point[1]) > dy) return kFALSE;
203  // then x
204  Double_t dx = 0.5*(fDx2*(point[2]+fDz)+fDx1*(fDz-point[2]))/fDz;
205  if (TMath::Abs(point[0]) > dx) return kFALSE;
206  return kTRUE;
207 }
208 
209 ////////////////////////////////////////////////////////////////////////////////
210 /// Compute distance from inside point to surface of the trd2
211 /// Boundary safe algorithm
212 
213 Double_t TGeoTrd2::DistFromInside(const Double_t *point, const Double_t *dir, Int_t iact, Double_t step, Double_t *safe) const
214 {
215  Double_t snxt = TGeoShape::Big();
216  if (iact<3 && safe) {
217  // compute safe distance
218  *safe = Safety(point, kTRUE);
219  if (iact==0) return TGeoShape::Big();
220  if (iact==1 && step<*safe) return TGeoShape::Big();
221  }
222 
223  Double_t fx = 0.5*(fDx1-fDx2)/fDz;
224  Double_t fy = 0.5*(fDy1-fDy2)/fDz;
225  Double_t cn;
226 
227  Double_t distx = 0.5*(fDx1+fDx2)-fx*point[2];
228  Double_t disty = 0.5*(fDy1+fDy2)-fy*point[2];
229 
230  //--- Compute distance to this shape
231  // first check if Z facettes are crossed
232  Double_t dist[3];
233  for (Int_t i=0; i<3; i++) dist[i]=TGeoShape::Big();
234  if (dir[2]<0) {
235  dist[0]=-(point[2]+fDz)/dir[2];
236  } else if (dir[2]>0) {
237  dist[0]=(fDz-point[2])/dir[2];
238  }
239  if (dist[0]<=0) return 0.0;
240  // now check X facettes
241  cn = -dir[0]+fx*dir[2];
242  if (cn>0) {
243  dist[1] = point[0]+distx;
244  if (dist[1]<=0) return 0.0;
245  dist[1] /= cn;
246  }
247  cn = dir[0]+fx*dir[2];
248  if (cn>0) {
249  Double_t s = distx-point[0];
250  if (s<=0) return 0.0;
251  s /= cn;
252  if (s<dist[1]) dist[1] = s;
253  }
254  // now check Y facettes
255  cn = -dir[1]+fy*dir[2];
256  if (cn>0) {
257  dist[2] = point[1]+disty;
258  if (dist[2]<=0) return 0.0;
259  dist[2] /= cn;
260  }
261  cn = dir[1]+fy*dir[2];
262  if (cn>0) {
263  Double_t s = disty-point[1];
264  if (s<=0) return 0.0;
265  s /= cn;
266  if (s<dist[2]) dist[2] = s;
267  }
268  snxt = dist[TMath::LocMin(3,dist)];
269  return snxt;
270 }
271 
272 ////////////////////////////////////////////////////////////////////////////////
273 /// Compute distance from outside point to surface of the trd2
274 /// Boundary safe algorithm
275 
276 Double_t TGeoTrd2::DistFromOutside(const Double_t *point, const Double_t *dir, Int_t iact, Double_t step, Double_t *safe) const
277 {
278  Double_t snxt = TGeoShape::Big();
279  if (iact<3 && safe) {
280  // compute safe distance
281  *safe = Safety(point, kFALSE);
282  if (iact==0) return TGeoShape::Big();
283  if (iact==1 && step<*safe) return TGeoShape::Big();
284  }
285  // find a visible face
286  Double_t xnew,ynew,znew;
287  Double_t fx = 0.5*(fDx1-fDx2)/fDz;
288  Double_t fy = 0.5*(fDy1-fDy2)/fDz;
289  Double_t cn;
290  // check visibility of X faces
291  Double_t distx = 0.5*(fDx1+fDx2)-fx*point[2];
292  Double_t disty = 0.5*(fDy1+fDy2)-fy*point[2];
293  Bool_t in = kTRUE;
294  Double_t safx = distx-TMath::Abs(point[0]);
295  Double_t safy = disty-TMath::Abs(point[1]);
296  Double_t safz = fDz-TMath::Abs(point[2]);
297  //--- Compute distance to this shape
298  // first check if Z facettes are crossed
299  if (point[2]<=-fDz) {
300  cn = -dir[2];
301  if (cn>=0) return TGeoShape::Big();
302  in = kFALSE;
303  snxt = (fDz+point[2])/cn;
304  // find extrapolated X and Y
305  xnew = point[0]+snxt*dir[0];
306  if (TMath::Abs(xnew) < fDx1) {
307  ynew = point[1]+snxt*dir[1];
308  if (TMath::Abs(ynew) < fDy1) return snxt;
309  }
310  } else if (point[2]>=fDz) {
311  cn = dir[2];
312  if (cn>=0) return TGeoShape::Big();
313  in = kFALSE;
314  snxt = (fDz-point[2])/cn;
315  // find extrapolated X and Y
316  xnew = point[0]+snxt*dir[0];
317  if (TMath::Abs(xnew) < fDx2) {
318  ynew = point[1]+snxt*dir[1];
319  if (TMath::Abs(ynew) < fDy2) return snxt;
320  }
321  }
322  // check if X facettes are crossed
323  if (point[0]<=-distx) {
324  cn = -dir[0]+fx*dir[2];
325  if (cn>=0) return TGeoShape::Big();
326  in = kFALSE;
327  snxt = (point[0]+distx)/cn;
328  // find extrapolated Y and Z
329  znew = point[2]+snxt*dir[2];
330  if (TMath::Abs(znew) < fDz) {
331  Double_t dy = 0.5*(fDy1+fDy2)-fy*znew;
332  ynew = point[1]+snxt*dir[1];
333  if (TMath::Abs(ynew) < dy) return snxt;
334  }
335  }
336  if (point[0]>=distx) {
337  cn = dir[0]+fx*dir[2];
338  if (cn>=0) return TGeoShape::Big();
339  in = kFALSE;
340  snxt = (distx-point[0])/cn;
341  // find extrapolated Y and Z
342  znew = point[2]+snxt*dir[2];
343  if (TMath::Abs(znew) < fDz) {
344  Double_t dy = 0.5*(fDy1+fDy2)-fy*znew;
345  ynew = point[1]+snxt*dir[1];
346  if (TMath::Abs(ynew) < dy) return snxt;
347  }
348  }
349  // finally check Y facettes
350  if (point[1]<=-disty) {
351  cn = -dir[1]+fy*dir[2];
352  in = kFALSE;
353  if (cn>=0) return TGeoShape::Big();
354  snxt = (point[1]+disty)/cn;
355  // find extrapolated X and Z
356  znew = point[2]+snxt*dir[2];
357  if (TMath::Abs(znew) < fDz) {
358  Double_t dx = 0.5*(fDx1+fDx2)-fx*znew;
359  xnew = point[0]+snxt*dir[0];
360  if (TMath::Abs(xnew) < dx) return snxt;
361  }
362  }
363  if (point[1]>=disty) {
364  cn = dir[1]+fy*dir[2];
365  if (cn>=0) return TGeoShape::Big();
366  in = kFALSE;
367  snxt = (disty-point[1])/cn;
368  // find extrapolated X and Z
369  znew = point[2]+snxt*dir[2];
370  if (TMath::Abs(znew) < fDz) {
371  Double_t dx = 0.5*(fDx1+fDx2)-fx*znew;
372  xnew = point[0]+snxt*dir[0];
373  if (TMath::Abs(xnew) < dx) return snxt;
374  }
375  }
376  if (!in) return TGeoShape::Big();
377  // Point actually inside
378  if (safz<safx && safz<safy) {
379  if (point[2]*dir[2]>=0) return TGeoShape::Big();
380  return 0.0;
381  }
382  if (safy<safx) {
383  cn = TMath::Sign(1.0,point[1])*dir[1]+fy*dir[2];
384  if (cn>=0) return TGeoShape::Big();
385  return 0.0;
386  }
387  cn = TMath::Sign(1.0,point[0])*dir[0]+fx*dir[2];
388  if (cn>=0) return TGeoShape::Big();
389  return 0.0;
390 }
391 
392 ////////////////////////////////////////////////////////////////////////////////
393 /// Get range of shape for a given axis.
394 
395 Double_t TGeoTrd2::GetAxisRange(Int_t iaxis, Double_t &xlo, Double_t &xhi) const
396 {
397  xlo = 0;
398  xhi = 0;
399  Double_t dx = 0;
400  switch (iaxis) {
401  case 3:
402  xlo = -fDz;
403  xhi = fDz;
404  dx = xhi-xlo;
405  return dx;
406  }
407  return dx;
408 }
409 
410 ////////////////////////////////////////////////////////////////////////////////
411 /// get the most visible corner from outside point and the normals
412 
413 void TGeoTrd2::GetVisibleCorner(const Double_t *point, Double_t *vertex, Double_t *normals) const
414 {
415  Double_t fx = 0.5*(fDx1-fDx2)/fDz;
416  Double_t fy = 0.5*(fDy1-fDy2)/fDz;
417  Double_t calf = 1./TMath::Sqrt(1.0+fx*fx);
418  Double_t salf = calf*fx;
419  Double_t cbet = 1./TMath::Sqrt(1.0+fy*fy);
420  Double_t sbet = cbet*fy;
421  // check visibility of X,Y faces
422  Double_t distx = fDx1-fx*(fDz+point[2]);
423  Double_t disty = fDy1-fy*(fDz+point[2]);
424  memset(normals, 0, 9*sizeof(Double_t));
425  TGeoTrd2 *trd2 = (TGeoTrd2*)this;
426  if (point[0]>distx) {
427  // hi x face visible
428  trd2->SetShapeBit(kGeoVisX);
429  normals[0]=calf;
430  normals[2]=salf;
431  } else {
432  trd2->SetShapeBit(kGeoVisX, kFALSE);
433  normals[0]=-calf;
434  normals[2]=salf;
435  }
436  if (point[1]>disty) {
437  // hi y face visible
438  trd2->SetShapeBit(kGeoVisY);
439  normals[4]=cbet;
440  normals[5]=sbet;
441  } else {
442  trd2->SetShapeBit(kGeoVisY, kFALSE);
443  normals[4]=-cbet;
444  normals[5]=sbet;
445  }
446  if (point[2]>fDz) {
447  // hi z face visible
448  trd2->SetShapeBit(kGeoVisZ);
449  normals[8]=1;
450  } else {
451  trd2->SetShapeBit(kGeoVisZ, kFALSE);
452  normals[8]=-1;
453  }
454  SetVertex(vertex);
455 }
456 
457 ////////////////////////////////////////////////////////////////////////////////
458 /// get the opposite corner of the intersected face
459 
460 void TGeoTrd2::GetOppositeCorner(const Double_t * /*point*/, Int_t inorm, Double_t *vertex, Double_t *normals) const
461 {
462  TGeoTrd2 *trd2 = (TGeoTrd2*)this;
463  if (inorm != 0) {
464  // change x face
465  trd2->SetShapeBit(kGeoVisX, !TestShapeBit(kGeoVisX));
466  normals[0]=-normals[0];
467  }
468  if (inorm != 1) {
469  // change y face
470  trd2->SetShapeBit(kGeoVisY, !TestShapeBit(kGeoVisY));
471  normals[4]=-normals[4];
472  }
473  if (inorm != 2) {
474  // hi z face visible
475  trd2->SetShapeBit(kGeoVisZ, !TestShapeBit(kGeoVisZ));
476  normals[8]=-normals[8];
477  }
478  SetVertex(vertex);
479 }
480 
481 ////////////////////////////////////////////////////////////////////////////////
482 /// Divide this trd2 shape belonging to volume "voldiv" into ndiv volumes
483 /// called divname, from start position with the given step. Only Z divisions
484 /// are supported. For Z divisions just return the pointer to the volume to be
485 /// divided. In case a wrong division axis is supplied, returns pointer to
486 /// volume that was divided.
487 
488 TGeoVolume *TGeoTrd2::Divide(TGeoVolume *voldiv, const char *divname, Int_t iaxis, Int_t ndiv,
489  Double_t start, Double_t step)
490 {
491  TGeoShape *shape; //--- shape to be created
492  TGeoVolume *vol; //--- division volume to be created
493  TGeoVolumeMulti *vmulti; //--- generic divided volume
494  TGeoPatternFinder *finder; //--- finder to be attached
495  TString opt = ""; //--- option to be attached
496  Double_t zmin, zmax, dx1n, dx2n, dy1n, dy2n;
497  Int_t id;
498  Double_t end = start+ndiv*step;
499  switch (iaxis) {
500  case 1:
501  Warning("Divide", "dividing a Trd2 on X not implemented");
502  return 0;
503  case 2:
504  Warning("Divide", "dividing a Trd2 on Y not implemented");
505  return 0;
506  case 3:
507  finder = new TGeoPatternZ(voldiv, ndiv, start, end);
508  vmulti = gGeoManager->MakeVolumeMulti(divname, voldiv->GetMedium());
509  voldiv->SetFinder(finder);
510  finder->SetDivIndex(voldiv->GetNdaughters());
511  for (id=0; id<ndiv; id++) {
512  zmin = start+id*step;
513  zmax = start+(id+1)*step;
514  dx1n = 0.5*(fDx1*(fDz-zmin)+fDx2*(fDz+zmin))/fDz;
515  dx2n = 0.5*(fDx1*(fDz-zmax)+fDx2*(fDz+zmax))/fDz;
516  dy1n = 0.5*(fDy1*(fDz-zmin)+fDy2*(fDz+zmin))/fDz;
517  dy2n = 0.5*(fDy1*(fDz-zmax)+fDy2*(fDz+zmax))/fDz;
518  shape = new TGeoTrd2(dx1n, dx2n, dy1n, dy2n, step/2.);
519  vol = new TGeoVolume(divname, shape, voldiv->GetMedium());
520  vmulti->AddVolume(vol);
521  opt = "Z";
522  voldiv->AddNodeOffset(vol, id, start+step/2+id*step, opt.Data());
523  ((TGeoNodeOffset*)voldiv->GetNodes()->At(voldiv->GetNdaughters()-1))->SetFinder(finder);
524  }
525  return vmulti;
526  default:
527  Error("Divide", "Wrong axis type for division");
528  return 0;
529  }
530 }
531 
532 ////////////////////////////////////////////////////////////////////////////////
533 /// Fill vector param[4] with the bounding cylinder parameters. The order
534 /// is the following : Rmin, Rmax, Phi1, Phi2
535 
536 void TGeoTrd2::GetBoundingCylinder(Double_t *param) const
537 {
538  TGeoBBox::GetBoundingCylinder(param);
539 }
540 
541 ////////////////////////////////////////////////////////////////////////////////
542 /// Fills real parameters of a positioned box inside this. Returns 0 if successful.
543 
544 Int_t TGeoTrd2::GetFittingBox(const TGeoBBox *parambox, TGeoMatrix *mat, Double_t &dx, Double_t &dy, Double_t &dz) const
545 {
546  dx=dy=dz=0;
547  if (mat->IsRotation()) {
548  Error("GetFittingBox", "cannot handle parametrized rotated volumes");
549  return 1; // ### rotation not accepted ###
550  }
551  //--> translate the origin of the parametrized box to the frame of this box.
552  Double_t origin[3];
553  mat->LocalToMaster(parambox->GetOrigin(), origin);
554  if (!Contains(origin)) {
555  Error("GetFittingBox", "wrong matrix - parametrized box is outside this");
556  return 1; // ### wrong matrix ###
557  }
558  //--> now we have to get the valid range for all parametrized axis
559  Double_t dd[3];
560  dd[0] = parambox->GetDX();
561  dd[1] = parambox->GetDY();
562  dd[2] = parambox->GetDZ();
563  //-> check if Z range is fixed
564  if (dd[2]<0) {
565  dd[2] = TMath::Min(origin[2]+fDz, fDz-origin[2]);
566  if (dd[2]<0) {
567  Error("GetFittingBox", "wrong matrix");
568  return 1;
569  }
570  }
571  if (dd[0]>=0 && dd[1]>=0) {
572  dx = dd[0];
573  dy = dd[1];
574  dz = dd[2];
575  return 0;
576  }
577  //-> check now range at Z = origin[2] +/- dd[2]
578  Double_t fx = 0.5*(fDx1-fDx2)/fDz;
579  Double_t fy = 0.5*(fDy1-fDy2)/fDz;
580  Double_t dx0 = 0.5*(fDx1+fDx2);
581  Double_t dy0 = 0.5*(fDy1+fDy2);
582  Double_t z=origin[2]-dd[2];
583  dd[0] = dx0-fx*z-origin[0];
584  dd[1] = dy0-fy*z-origin[1];
585  z=origin[2]+dd[2];
586  dd[0] = TMath::Min(dd[0], dx0-fx*z-origin[0]);
587  dd[1] = TMath::Min(dd[1], dy0-fy*z-origin[1]);
588  if (dd[0]<0 || dd[1]<0) {
589  Error("GetFittingBox", "wrong matrix");
590  return 1;
591  }
592  dx = dd[0];
593  dy = dd[1];
594  dz = dd[2];
595  return 0;
596 }
597 
598 ////////////////////////////////////////////////////////////////////////////////
599 /// in case shape has some negative parameters, these has to be computed
600 /// in order to fit the mother
601 
602 TGeoShape *TGeoTrd2::GetMakeRuntimeShape(TGeoShape *mother, TGeoMatrix * /*mat*/) const
603 {
604  if (!TestShapeBit(kGeoRunTimeShape)) return 0;
605  if (!mother->TestShapeBit(kGeoTrd2)) {
606  Error("GetMakeRuntimeShape", "invalid mother");
607  return 0;
608  }
609  Double_t dx1, dx2, dy1, dy2, dz;
610  if (fDx1<0) dx1=((TGeoTrd2*)mother)->GetDx1();
611  else dx1=fDx1;
612  if (fDx2<0) dx2=((TGeoTrd2*)mother)->GetDx2();
613  else dx2=fDx2;
614  if (fDy1<0) dy1=((TGeoTrd2*)mother)->GetDy1();
615  else dy1=fDy1;
616  if (fDy2<0) dy2=((TGeoTrd2*)mother)->GetDy2();
617  else dy2=fDy2;
618  if (fDz<0) dz=((TGeoTrd2*)mother)->GetDz();
619  else dz=fDz;
620 
621  return (new TGeoTrd2(dx1, dx2, dy1, dy2, dz));
622 }
623 
624 ////////////////////////////////////////////////////////////////////////////////
625 /// print shape parameters
626 
627 void TGeoTrd2::InspectShape() const
628 {
629  printf("*** Shape %s: TGeoTrd2 ***\n", GetName());
630  printf(" dx1 = %11.5f\n", fDx1);
631  printf(" dx2 = %11.5f\n", fDx2);
632  printf(" dy1 = %11.5f\n", fDy1);
633  printf(" dy2 = %11.5f\n", fDy2);
634  printf(" dz = %11.5f\n", fDz);
635  printf(" Bounding box:\n");
636  TGeoBBox::InspectShape();
637 }
638 
639 ////////////////////////////////////////////////////////////////////////////////
640 /// computes the closest distance from given point to this shape, according
641 /// to option. The matching point on the shape is stored in spoint.
642 
643 Double_t TGeoTrd2::Safety(const Double_t *point, Bool_t in) const
644 {
645  Double_t saf[3];
646  //--- Compute safety first
647  // check Z facettes
648  saf[0] = fDz-TMath::Abs(point[2]);
649  Double_t fx = 0.5*(fDx1-fDx2)/fDz;
650  Double_t calf = 1./TMath::Sqrt(1.0+fx*fx);
651  // check X facettes
652  Double_t distx = 0.5*(fDx1+fDx2)-fx*point[2];
653  if (distx<0) saf[1]=TGeoShape::Big();
654  else saf[1]=(distx-TMath::Abs(point[0]))*calf;
655 
656  Double_t fy = 0.5*(fDy1-fDy2)/fDz;
657  calf = 1./TMath::Sqrt(1.0+fy*fy);
658  // check Y facettes
659  distx = 0.5*(fDy1+fDy2)-fy*point[2];
660  if (distx<0) saf[2]=TGeoShape::Big();
661  else saf[2]=(distx-TMath::Abs(point[1]))*calf;
662 
663  if (in) return saf[TMath::LocMin(3,saf)];
664  for (Int_t i=0; i<3; i++) saf[i]=-saf[i];
665  return saf[TMath::LocMax(3,saf)];
666 }
667 
668 ////////////////////////////////////////////////////////////////////////////////
669 /// Save a primitive as a C++ statement(s) on output stream "out".
670 
671 void TGeoTrd2::SavePrimitive(std::ostream &out, Option_t * /*option*/ /*= ""*/)
672 {
673  if (TObject::TestBit(kGeoSavePrimitive)) return;
674  out << " // Shape: " << GetName() << " type: " << ClassName() << std::endl;
675  out << " dx1 = " << fDx1 << ";" << std::endl;
676  out << " dx2 = " << fDx2 << ";" << std::endl;
677  out << " dy1 = " << fDy1 << ";" << std::endl;
678  out << " dy2 = " << fDy2 << ";" << std::endl;
679  out << " dz = " << fDZ << ";" << std::endl;
680  out << " TGeoShape *" << GetPointerName() << " = new TGeoTrd2(\"" << GetName() << "\", dx1,dx2,dy1,dy2,dz);" << std::endl;
681  TObject::SetBit(TGeoShape::kGeoSavePrimitive);
682 }
683 
684 ////////////////////////////////////////////////////////////////////////////////
685 /// set arb8 params in one step :
686 
687 void TGeoTrd2::SetDimensions(Double_t *param)
688 {
689  fDx1 = param[0];
690  fDx2 = param[1];
691  fDy1 = param[2];
692  fDy2 = param[3];
693  fDz = param[4];
694  ComputeBBox();
695 }
696 
697 ////////////////////////////////////////////////////////////////////////////////
698 /// create trd2 mesh points
699 
700 void TGeoTrd2::SetPoints(Double_t *points) const
701 {
702  if (!points) return;
703  points[ 0] = -fDx1; points[ 1] = -fDy1; points[ 2] = -fDz;
704  points[ 3] = -fDx1; points[ 4] = fDy1; points[ 5] = -fDz;
705  points[ 6] = fDx1; points[ 7] = fDy1; points[ 8] = -fDz;
706  points[ 9] = fDx1; points[10] = -fDy1; points[11] = -fDz;
707  points[12] = -fDx2; points[13] = -fDy2; points[14] = fDz;
708  points[15] = -fDx2; points[16] = fDy2; points[17] = fDz;
709  points[18] = fDx2; points[19] = fDy2; points[20] = fDz;
710  points[21] = fDx2; points[22] = -fDy2; points[23] = fDz;
711 }
712 
713 ////////////////////////////////////////////////////////////////////////////////
714 /// create trd2 mesh points
715 
716 void TGeoTrd2::SetPoints(Float_t *points) const
717 {
718  if (!points) return;
719  points[ 0] = -fDx1; points[ 1] = -fDy1; points[ 2] = -fDz;
720  points[ 3] = -fDx1; points[ 4] = fDy1; points[ 5] = -fDz;
721  points[ 6] = fDx1; points[ 7] = fDy1; points[ 8] = -fDz;
722  points[ 9] = fDx1; points[10] = -fDy1; points[11] = -fDz;
723  points[12] = -fDx2; points[13] = -fDy2; points[14] = fDz;
724  points[15] = -fDx2; points[16] = fDy2; points[17] = fDz;
725  points[18] = fDx2; points[19] = fDy2; points[20] = fDz;
726  points[21] = fDx2; points[22] = -fDy2; points[23] = fDz;
727 }
728 
729 ////////////////////////////////////////////////////////////////////////////////
730 /// set vertex of a corner according to visibility flags
731 
732 void TGeoTrd2::SetVertex(Double_t *vertex) const
733 {
734  if (TestShapeBit(kGeoVisX)) {
735  if (TestShapeBit(kGeoVisZ)) {
736  vertex[0] = fDx2;
737  vertex[2] = fDz;
738  vertex[1] = (TestShapeBit(kGeoVisY))?fDy2:-fDy2;
739  } else {
740  vertex[0] = fDx1;
741  vertex[2] = -fDz;
742  vertex[1] = (TestShapeBit(kGeoVisY))?fDy1:-fDy1;
743  }
744  } else {
745  if (TestShapeBit(kGeoVisZ)) {
746  vertex[0] = -fDx2;
747  vertex[2] = fDz;
748  vertex[1] = (TestShapeBit(kGeoVisY))?fDy2:-fDy2;
749  } else {
750  vertex[0] = -fDx1;
751  vertex[2] = -fDz;
752  vertex[1] = (TestShapeBit(kGeoVisY))?fDy1:-fDy1;
753  }
754  }
755 }
756 
757 ////////////////////////////////////////////////////////////////////////////////
758 /// fill size of this 3-D object
759 
760 void TGeoTrd2::Sizeof3D() const
761 {
762  TGeoBBox::Sizeof3D();
763 }
764 
765 ////////////////////////////////////////////////////////////////////////////////
766 /// Check the inside status for each of the points in the array.
767 /// Input: Array of point coordinates + vector size
768 /// Output: Array of Booleans for the inside of each point
769 
770 void TGeoTrd2::Contains_v(const Double_t *points, Bool_t *inside, Int_t vecsize) const
771 {
772  for (Int_t i=0; i<vecsize; i++) inside[i] = Contains(&points[3*i]);
773 }
774 
775 ////////////////////////////////////////////////////////////////////////////////
776 /// Compute the normal for an array o points so that norm.dot.dir is positive
777 /// Input: Arrays of point coordinates and directions + vector size
778 /// Output: Array of normal directions
779 
780 void TGeoTrd2::ComputeNormal_v(const Double_t *points, const Double_t *dirs, Double_t *norms, Int_t vecsize)
781 {
782  for (Int_t i=0; i<vecsize; i++) ComputeNormal(&points[3*i], &dirs[3*i], &norms[3*i]);
783 }
784 
785 ////////////////////////////////////////////////////////////////////////////////
786 /// Compute distance from array of input points having directions specified by dirs. Store output in dists
787 
788 void TGeoTrd2::DistFromInside_v(const Double_t *points, const Double_t *dirs, Double_t *dists, Int_t vecsize, Double_t* step) const
789 {
790  for (Int_t i=0; i<vecsize; i++) dists[i] = DistFromInside(&points[3*i], &dirs[3*i], 3, step[i]);
791 }
792 
793 ////////////////////////////////////////////////////////////////////////////////
794 /// Compute distance from array of input points having directions specified by dirs. Store output in dists
795 
796 void TGeoTrd2::DistFromOutside_v(const Double_t *points, const Double_t *dirs, Double_t *dists, Int_t vecsize, Double_t* step) const
797 {
798  for (Int_t i=0; i<vecsize; i++) dists[i] = DistFromOutside(&points[3*i], &dirs[3*i], 3, step[i]);
799 }
800 
801 ////////////////////////////////////////////////////////////////////////////////
802 /// Compute safe distance from each of the points in the input array.
803 /// Input: Array of point coordinates, array of statuses for these points, size of the arrays
804 /// Output: Safety values
805 
806 void TGeoTrd2::Safety_v(const Double_t *points, const Bool_t *inside, Double_t *safe, Int_t vecsize) const
807 {
808  for (Int_t i=0; i<vecsize; i++) safe[i] = Safety(&points[3*i], inside[i]);
809 }