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