Logo ROOT  
Reference Guide
TGeoVoxelFinder.cxx
Go to the documentation of this file.
1 // @(#)root/geom:$Id$
2 // Author: Andrei Gheata 04/02/02
3 
4 /*************************************************************************
5  * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers. *
6  * All rights reserved. *
7  * *
8  * For the licensing terms see $ROOTSYS/LICENSE. *
9  * For the list of contributors see $ROOTSYS/README/CREDITS. *
10  *************************************************************************/
11 
12 /** \class TGeoVoxelFinder
13 \ingroup Geometry_classes
14 
15 Finder class handling voxels.
16 
17 Full description with examples and pictures
18 
19 \image html geom_t_finder.png
20 \image html geom_t_voxelfind.png
21 \image html geom_t_voxtree.png
22 */
23 
24 #include "TGeoVoxelFinder.h"
25 
26 #include "TBuffer.h"
27 #include "TMath.h"
28 #include "TGeoMatrix.h"
29 #include "TGeoBBox.h"
30 #include "TGeoNode.h"
31 #include "TGeoManager.h"
32 #include "TGeoStateInfo.h"
33 
35 
36 ////////////////////////////////////////////////////////////////////////////////
37 /// Default constructor
38 
40 {
41  fVolume = 0;
42  fNboxes = 0;
43  fIbx = 0;
44  fIby = 0;
45  fIbz = 0;
46  fNox = 0;
47  fNoy = 0;
48  fNoz = 0;
49  fNex = 0;
50  fNey = 0;
51  fNez = 0;
52  fNx = 0;
53  fNy = 0;
54  fNz = 0;
55  fBoxes = 0;
56  fXb = 0;
57  fYb = 0;
58  fZb = 0;
59  fOBx = 0;
60  fOBy = 0;
61  fOBz = 0;
62  fOEx = 0;
63  fOEy = 0;
64  fOEz = 0;
65  fIndcX = 0;
66  fIndcY = 0;
67  fIndcZ = 0;
68  fExtraX = 0;
69  fExtraY = 0;
70  fExtraZ = 0;
71  fNsliceX = 0;
72  fNsliceY = 0;
73  fNsliceZ = 0;
74  memset(fPriority, 0, 3*sizeof(Int_t));
76 }
77 ////////////////////////////////////////////////////////////////////////////////
78 /// Default constructor
79 
81 {
82  if (!vol) {
83  Fatal("TGeoVoxelFinder", "Null pointer for volume");
84  return; // To make code checkers happy
85  }
86  fVolume = vol;
88  fNboxes = 0;
89  fIbx = 0;
90  fIby = 0;
91  fIbz = 0;
92  fNox = 0;
93  fNoy = 0;
94  fNoz = 0;
95  fNex = 0;
96  fNey = 0;
97  fNez = 0;
98  fNx = 0;
99  fNy = 0;
100  fNz = 0;
101  fBoxes = 0;
102  fXb = 0;
103  fYb = 0;
104  fZb = 0;
105  fOBx = 0;
106  fOBy = 0;
107  fOBz = 0;
108  fOEx = 0;
109  fOEy = 0;
110  fOEz = 0;
111  fIndcX = 0;
112  fIndcY = 0;
113  fIndcZ = 0;
114  fExtraX = 0;
115  fExtraY = 0;
116  fExtraZ = 0;
117  fNsliceX = 0;
118  fNsliceY = 0;
119  fNsliceZ = 0;
120  memset(fPriority, 0, 3*sizeof(Int_t));
121  SetNeedRebuild();
122 }
123 
124 ////////////////////////////////////////////////////////////////////////////////
125 /// Destructor
126 
128 {
129  if (fOBx) delete [] fOBx;
130  if (fOBy) delete [] fOBy;
131  if (fOBz) delete [] fOBz;
132  if (fOEx) delete [] fOEx;
133  if (fOEy) delete [] fOEy;
134  if (fOEz) delete [] fOEz;
135 // printf("OBx OBy OBz...\n");
136  if (fBoxes) delete [] fBoxes;
137 // printf("boxes...\n");
138  if (fXb) delete [] fXb;
139  if (fYb) delete [] fYb;
140  if (fZb) delete [] fZb;
141 // printf("Xb Yb Zb...\n");
142  if (fNsliceX) delete [] fNsliceX;
143  if (fNsliceY) delete [] fNsliceY;
144  if (fNsliceZ) delete [] fNsliceZ;
145  if (fIndcX) delete [] fIndcX;
146  if (fIndcY) delete [] fIndcY;
147  if (fIndcZ) delete [] fIndcZ;
148  if (fExtraX) delete [] fExtraX;
149  if (fExtraY) delete [] fExtraY;
150  if (fExtraZ) delete [] fExtraZ;
151 // printf("IndX IndY IndZ...\n");
152 }
153 
154 ////////////////////////////////////////////////////////////////////////////////
155 
157 {
158  return td.fVoxNcandidates;
159 }
160 
161 ////////////////////////////////////////////////////////////////////////////////
162 
164 {
165  nelem = td.fVoxNcandidates;
166  return td.fVoxCheckList;
167 }
168 
169 ////////////////////////////////////////////////////////////////////////////////
170 /// build the array of bounding boxes of the nodes inside
171 
173 {
174  Int_t nd = fVolume->GetNdaughters();
175  if (!nd) return;
176  Int_t id;
177  TGeoNode *node;
178  if (fBoxes) delete [] fBoxes;
179  fNboxes = 6*nd;
180  fBoxes = new Double_t[fNboxes];
181  Double_t vert[24] = {0};
182  Double_t pt[3] = {0};
183  Double_t xyz[6] = {0};
184 // printf("boundaries for %s :\n", GetName());
185  TGeoBBox *box = 0;
186  for (id=0; id<nd; id++) {
187  node = fVolume->GetNode(id);
188 // if (!strcmp(node->ClassName(), "TGeoNodeOffset") continue;
189  box = (TGeoBBox*)node->GetVolume()->GetShape();
190  box->SetBoxPoints(&vert[0]);
191  for (Int_t point=0; point<8; point++) {
192  DaughterToMother(id, &vert[3*point], &pt[0]);
193  if (!point) {
194  xyz[0] = xyz[1] = pt[0];
195  xyz[2] = xyz[3] = pt[1];
196  xyz[4] = xyz[5] = pt[2];
197  continue;
198  }
199  for (Int_t j=0; j<3; j++) {
200  if (pt[j] < xyz[2*j]) xyz[2*j]=pt[j];
201  if (pt[j] > xyz[2*j+1]) xyz[2*j+1]=pt[j];
202  }
203  }
204  fBoxes[6*id] = 0.5*(xyz[1]-xyz[0]); // dX
205  fBoxes[6*id+1] = 0.5*(xyz[3]-xyz[2]); // dY
206  fBoxes[6*id+2] = 0.5*(xyz[5]-xyz[4]); // dZ
207  fBoxes[6*id+3] = 0.5*(xyz[0]+xyz[1]); // Ox
208  fBoxes[6*id+4] = 0.5*(xyz[2]+xyz[3]); // Oy
209  fBoxes[6*id+5] = 0.5*(xyz[4]+xyz[5]); // Oz
210  }
211 }
212 
213 ////////////////////////////////////////////////////////////////////////////////
214 /// convert a point from the local reference system of node id to reference
215 /// system of mother volume
216 
217 void TGeoVoxelFinder::DaughterToMother(Int_t id, const Double_t *local, Double_t *master) const
218 {
219  TGeoMatrix *mat = fVolume->GetNode(id)->GetMatrix();
220  if (!mat) memcpy(master,local,3*sizeof(Double_t));
221  else mat->LocalToMaster(local, master);
222 }
223 ////////////////////////////////////////////////////////////////////////////////
224 /// Computes squared distance from POINT to the voxel(s) containing node INODE. Returns 0
225 /// if POINT inside voxel(s).
226 
227 Bool_t TGeoVoxelFinder::IsSafeVoxel(const Double_t *point, Int_t inode, Double_t minsafe) const
228 {
229  if (NeedRebuild()) {
230  TGeoVoxelFinder *vox = (TGeoVoxelFinder*)this;
231  vox->Voxelize();
233  }
234  Double_t dxyz, minsafe2=minsafe*minsafe;
235  Int_t ist = 6*inode;
236  Int_t i;
237  Double_t rsq = 0;
238  for (i=0; i<3; i++) {
239  dxyz = TMath::Abs(point[i]-fBoxes[ist+i+3])-fBoxes[ist+i];
240  if (dxyz>-1E-6) rsq+=dxyz*dxyz;
241  if (rsq > minsafe2*(1.+TGeoShape::Tolerance())) return kTRUE;
242  }
243  return kFALSE;
244 }
245 
246 ////////////////////////////////////////////////////////////////////////////////
247 /// Compute voxelization efficiency.
248 
250 {
251  printf("Voxelization efficiency for %s\n", fVolume->GetName());
252  if (NeedRebuild()) {
253  Voxelize();
255  }
257  Double_t eff = 0;
258  Double_t effslice = 0;
259  Int_t id;
260  if (fPriority[0]) {
261  for (id=0; id<fIbx-1; id++) { // loop on boundaries
262  effslice += fNsliceX[id];
263  }
264  if (!TGeoShape::IsSameWithinTolerance(effslice,0)) effslice = nd/effslice;
265  else printf("Woops : slice X\n");
266  }
267  printf("X efficiency : %g\n", effslice);
268  eff += effslice;
269  effslice = 0;
270  if (fPriority[1]) {
271  for (id=0; id<fIby-1; id++) { // loop on boundaries
272  effslice += fNsliceY[id];
273  }
274  if (!TGeoShape::IsSameWithinTolerance(effslice,0)) effslice = nd/effslice;
275  else printf("Woops : slice X\n");
276  }
277  printf("Y efficiency : %g\n", effslice);
278  eff += effslice;
279  effslice = 0;
280  if (fPriority[2]) {
281  for (id=0; id<fIbz-1; id++) { // loop on boundaries
282  effslice += fNsliceZ[id];
283  }
284  if (!TGeoShape::IsSameWithinTolerance(effslice,0)) effslice = nd/effslice;
285  else printf("Woops : slice X\n");
286  }
287  printf("Z efficiency : %g\n", effslice);
288  eff += effslice;
289  eff /= 3.;
290  printf("Total efficiency : %g\n", eff);
291  return eff;
292 }
293 ////////////////////////////////////////////////////////////////////////////////
294 /// create the list of nodes for which the bboxes overlap with inode's bbox
295 
297 {
298  if (!fBoxes) return;
299  Double_t xmin, xmax, ymin, ymax, zmin, zmax;
300  Double_t xmin1, xmax1, ymin1, ymax1, zmin1, zmax1;
301  Double_t ddx1, ddx2;
302  Int_t nd = fVolume->GetNdaughters();
303  Int_t *ovlps = 0;
304  Int_t *otmp = new Int_t[nd-1];
305  Int_t novlp = 0;
306  TGeoNode *node = fVolume->GetNode(inode);
307  xmin = fBoxes[6*inode+3] - fBoxes[6*inode];
308  xmax = fBoxes[6*inode+3] + fBoxes[6*inode];
309  ymin = fBoxes[6*inode+4] - fBoxes[6*inode+1];
310  ymax = fBoxes[6*inode+4] + fBoxes[6*inode+1];
311  zmin = fBoxes[6*inode+5] - fBoxes[6*inode+2];
312  zmax = fBoxes[6*inode+5] + fBoxes[6*inode+2];
313  // loop on brothers
314  for (Int_t ib=0; ib<nd; ib++) {
315  if (ib == inode) continue; // everyone overlaps with itself
316  xmin1 = fBoxes[6*ib+3] - fBoxes[6*ib];
317  xmax1 = fBoxes[6*ib+3] + fBoxes[6*ib];
318  ymin1 = fBoxes[6*ib+4] - fBoxes[6*ib+1];
319  ymax1 = fBoxes[6*ib+4] + fBoxes[6*ib+1];
320  zmin1 = fBoxes[6*ib+5] - fBoxes[6*ib+2];
321  zmax1 = fBoxes[6*ib+5] + fBoxes[6*ib+2];
322 
323  ddx1 = xmax-xmin1;
324  ddx2 = xmax1-xmin;
325  if (ddx1*ddx2 <= 0.) continue;
326  ddx1 = ymax-ymin1;
327  ddx2 = ymax1-ymin;
328  if (ddx1*ddx2 <= 0.) continue;
329  ddx1 = zmax-zmin1;
330  ddx2 = zmax1-zmin;
331  if (ddx1*ddx2 <= 0.) continue;
332  otmp[novlp++] = ib;
333  }
334  if (!novlp) {
335  delete [] otmp;
336  node->SetOverlaps(ovlps, 0);
337  return;
338  }
339  ovlps = new Int_t[novlp];
340  memcpy(ovlps, otmp, novlp*sizeof(Int_t));
341  delete [] otmp;
342  node->SetOverlaps(ovlps, novlp);
343 }
344 
345 ////////////////////////////////////////////////////////////////////////////////
346 /// Get indices for current slices on x, y, z
347 
349 {
350  td.fVoxSlices[0] = -2; // -2 means 'all daughters in slice'
351  td.fVoxSlices[1] = -2;
352  td.fVoxSlices[2] = -2;
353  Bool_t flag=kTRUE;
354  if (fPriority[0]) {
355  td.fVoxSlices[0] = TMath::BinarySearch(fIbx, fXb, point[0]);
356  if ((td.fVoxSlices[0]<0) || (td.fVoxSlices[0]>=fIbx-1)) {
357  // outside slices
358  flag=kFALSE;
359  } else {
360  if (fPriority[0]==2) {
361  // nothing in current slice
362  if (!fNsliceX[td.fVoxSlices[0]]) flag = kFALSE;
363  }
364  }
365  }
366  if (fPriority[1]) {
367  td.fVoxSlices[1] = TMath::BinarySearch(fIby, fYb, point[1]);
368  if ((td.fVoxSlices[1]<0) || (td.fVoxSlices[1]>=fIby-1)) {
369  // outside slices
370  flag=kFALSE;
371  } else {
372  if (fPriority[1]==2) {
373  // nothing in current slice
374  if (!fNsliceY[td.fVoxSlices[1]]) flag = kFALSE;
375  }
376  }
377  }
378  if (fPriority[2]) {
379  td.fVoxSlices[2] = TMath::BinarySearch(fIbz, fZb, point[2]);
380  if ((td.fVoxSlices[2]<0) || (td.fVoxSlices[2]>=fIbz-1)) return kFALSE;
381  if (fPriority[2]==2) {
382  // nothing in current slice
383  if (!fNsliceZ[td.fVoxSlices[2]]) return kFALSE;
384  }
385  }
386  return flag;
387 }
388 
389 ////////////////////////////////////////////////////////////////////////////////
390 /// Return the list of extra candidates in a given X slice compared to
391 /// another (left or right)
392 
393 Int_t *TGeoVoxelFinder::GetExtraX(Int_t islice, Bool_t left, Int_t &nextra) const
394 {
395  Int_t *list = 0;
396  nextra = 0;
397  if (fPriority[0]!=2) return list;
398  if (left) {
399  nextra = fExtraX[fOEx[islice]];
400  list = &fExtraX[fOEx[islice]+2];
401  } else {
402  nextra = fExtraX[fOEx[islice]+1];
403  list = &fExtraX[fOEx[islice]+2+fExtraX[fOEx[islice]]];
404  }
405  return list;
406 }
407 
408 ////////////////////////////////////////////////////////////////////////////////
409 /// Return the list of extra candidates in a given Y slice compared to
410 /// another (left or right)
411 
412 Int_t *TGeoVoxelFinder::GetExtraY(Int_t islice, Bool_t left, Int_t &nextra) const
413 {
414  Int_t *list = 0;
415  nextra = 0;
416  if (fPriority[1]!=2) return list;
417  if (left) {
418  nextra = fExtraY[fOEy[islice]];
419  list = &fExtraY[fOEy[islice]+2];
420  } else {
421  nextra = fExtraY[fOEy[islice]+1];
422  list = &fExtraY[fOEy[islice]+2+fExtraY[fOEy[islice]]];
423  }
424  return list;
425 }
426 
427 ////////////////////////////////////////////////////////////////////////////////
428 /// Return the list of extra candidates in a given Z slice compared to
429 /// another (left or right)
430 
431 Int_t *TGeoVoxelFinder::GetExtraZ(Int_t islice, Bool_t left, Int_t &nextra) const
432 {
433  Int_t *list = 0;
434  nextra = 0;
435  if (fPriority[2]!=2) return list;
436  if (left) {
437  nextra = fExtraZ[fOEz[islice]];
438  list = &fExtraZ[fOEz[islice]+2];
439  } else {
440  nextra = fExtraZ[fOEz[islice]+1];
441  list = &fExtraZ[fOEz[islice]+2+fExtraZ[fOEz[islice]]];
442  }
443  return list;
444 }
445 
446 ////////////////////////////////////////////////////////////////////////////////
447 /// Get extra candidates that are not contained in current check list
448 
450 {
451  td.fVoxNcandidates = 0;
452  Int_t icand;
453  UInt_t bitnumber, loc;
454  UChar_t bit, byte;
455  for (icand=0; icand<ncheck; icand++) {
456  bitnumber = (UInt_t)list[icand];
457  loc = bitnumber>>3;
458  bit = bitnumber%8;
459  byte = (~td.fVoxBits1[loc]) & (1<<bit);
460  if (byte) td.fVoxCheckList[td.fVoxNcandidates++]=list[icand];
461  }
462  ncheck = td.fVoxNcandidates;
463  return td.fVoxCheckList;
464 }
465 
466 ////////////////////////////////////////////////////////////////////////////////
467 /// Get extra candidates that are contained in array1 but not in current check list
468 
470 {
471  td.fVoxNcandidates = 0;
472  Int_t icand;
473  UInt_t bitnumber, loc;
474  UChar_t bit, byte;
475  for (icand=0; icand<ncheck; icand++) {
476  bitnumber = (UInt_t)list[icand];
477  loc = bitnumber>>3;
478  bit = bitnumber%8;
479  byte = (~td.fVoxBits1[loc]) & array1[loc] & (1<<bit);
480  if (byte) td.fVoxCheckList[td.fVoxNcandidates++]=list[icand];
481  }
482  ncheck = td.fVoxNcandidates;
483  return td.fVoxCheckList;
484 }
485 
486 ////////////////////////////////////////////////////////////////////////////////
487 /// Get extra candidates that are contained in array1 but not in current check list
488 
489 Int_t *TGeoVoxelFinder::GetValidExtra(Int_t /*n1*/, UChar_t *array1, Int_t /*n2*/, UChar_t *array2, Int_t *list, Int_t &ncheck, TGeoStateInfo &td)
490 {
491  td.fVoxNcandidates = 0;
492  Int_t icand;
493  UInt_t bitnumber, loc;
494  UChar_t bit, byte;
495  for (icand=0; icand<ncheck; icand++) {
496  bitnumber = (UInt_t)list[icand];
497  loc = bitnumber>>3;
498  bit = bitnumber%8;
499  byte = (~td.fVoxBits1[loc]) & array1[loc] & array2[loc] & (1<<bit);
500  if (byte) td.fVoxCheckList[td.fVoxNcandidates++]=list[icand];
501  }
502  ncheck = td.fVoxNcandidates;
503  return td.fVoxCheckList;
504 }
505 
506 ////////////////////////////////////////////////////////////////////////////////
507 /// Returns list of new candidates in next voxel. If NULL, nowhere to
508 /// go next.
509 
511 {
512  if (NeedRebuild()) {
513  Voxelize();
515  }
516  ncheck = 0;
517  if (td.fVoxLimits[0]<0) return 0;
518  if (td.fVoxLimits[1]<0) return 0;
519  if (td.fVoxLimits[2]<0) return 0;
520  Int_t dind[3]; // new slices
521  //---> start from old slices
522  memcpy(&dind[0], &td.fVoxSlices[0], 3*sizeof(Int_t));
523  Double_t dmin[3]; // distances to get to next X,Y, Z slices.
524  dmin[0] = dmin[1] = dmin[2] = TGeoShape::Big();
525  //---> max. possible step to be considered
527 // printf("1- maxstep=%g\n", maxstep);
528  Bool_t isXlimit=kFALSE, isYlimit=kFALSE, isZlimit=kFALSE;
529  Bool_t isForcedX=kFALSE, isForcedY=kFALSE, isForcedZ=kFALSE;
530  Double_t dforced[3];
531  dforced[0] = dforced[1] = dforced[2] = TGeoShape::Big();
532  Int_t iforced = 0;
533  //
534  //---> work on X
535  if (fPriority[0] && td.fVoxInc[0]) {
536  //---> increment/decrement slice
537  dind[0] += td.fVoxInc[0];
538  if (td.fVoxInc[0]==1) {
539  if (dind[0]<0 || dind[0]>fIbx-1) return 0; // outside range
540  dmin[0] = (fXb[dind[0]]-point[0])*td.fVoxInvdir[0];
541  } else {
542  if (td.fVoxSlices[0]<0 || td.fVoxSlices[0]>fIbx-1) return 0; // outside range
543  dmin[0] = (fXb[td.fVoxSlices[0]]-point[0])*td.fVoxInvdir[0];
544  }
545  isXlimit = (dmin[0]>maxstep)?kTRUE:kFALSE;
546 // printf("---> X : priority=%i, slice=%i/%i inc=%i\n",
547 // fPriority[0], td.fVoxSlices[0], fIbx-2, td.fVoxInc[0]);
548 // printf("2- step to next X (%i) = %g\n", (Int_t)isXlimit, dmin[0]);
549  //---> check if propagation to next slice on this axis is forced
550  if ((td.fVoxSlices[0]==-1) || (td.fVoxSlices[0]==fIbx-1)) {
551  isForcedX = kTRUE;
552  dforced[0] = dmin[0];
553  iforced++;
554 // printf(" FORCED 1\n");
555  if (isXlimit) return 0;
556  } else {
557  if (fPriority[0]==2) {
558  // if no candidates in current slice, force next slice
559  if (fNsliceX[td.fVoxSlices[0]]==0) {
560  isForcedX = kTRUE;
561  dforced[0] = dmin[0];
562  iforced++;
563 // printf(" FORCED 2\n");
564  if (isXlimit) return 0;
565  }
566  }
567  }
568  } else {
569  // no slices on this axis -> bounding box limit
570 // printf(" No slice on X\n");
571  dmin[0] = td.fVoxLimits[0];
572  isXlimit = kTRUE;
573  }
574  //---> work on Y
575  if (fPriority[1] && td.fVoxInc[1]) {
576  //---> increment/decrement slice
577  dind[1] += td.fVoxInc[1];
578  if (td.fVoxInc[1]==1) {
579  if (dind[1]<0 || dind[1]>fIby-1) return 0; // outside range
580  dmin[1] = (fYb[dind[1]]-point[1])*td.fVoxInvdir[1];
581  } else {
582  if (td.fVoxSlices[1]<0 || td.fVoxSlices[1]>fIby-1) return 0; // outside range
583  dmin[1] = (fYb[td.fVoxSlices[1]]-point[1])*td.fVoxInvdir[1];
584  }
585  isYlimit = (dmin[1]>maxstep)?kTRUE:kFALSE;
586 // printf("---> Y : priority=%i, slice=%i/%i inc=%i\n",
587 // fPriority[1], td.fVoxSlices[1], fIby-2, td.fVoxInc[1]);
588 // printf("3- step to next Y (%i) = %g\n", (Int_t)isYlimit, dmin[1]);
589 
590  //---> check if propagation to next slice on this axis is forced
591  if ((td.fVoxSlices[1]==-1) || (td.fVoxSlices[1]==fIby-1)) {
592  isForcedY = kTRUE;
593  dforced[1] = dmin[1];
594  iforced++;
595 // printf(" FORCED 1\n");
596  if (isYlimit) return 0;
597  } else {
598  if (fPriority[1]==2) {
599  // if no candidates in current slice, force next slice
600  if (fNsliceY[td.fVoxSlices[1]]==0) {
601  isForcedY = kTRUE;
602  dforced[1] = dmin[1];
603  iforced++;
604 // printf(" FORCED 2\n");
605  if (isYlimit) return 0;
606  }
607  }
608  }
609  } else {
610  // no slices on this axis -> bounding box limit
611 // printf(" No slice on Y\n");
612  dmin[1] = td.fVoxLimits[1];
613  isYlimit = kTRUE;
614  }
615  //---> work on Z
616  if (fPriority[2] && td.fVoxInc[2]) {
617  //---> increment/decrement slice
618  dind[2] += td.fVoxInc[2];
619  if (td.fVoxInc[2]==1) {
620  if (dind[2]<0 || dind[2]>fIbz-1) return 0; // outside range
621  dmin[2] = (fZb[dind[2]]-point[2])*td.fVoxInvdir[2];
622  } else {
623  if (td.fVoxSlices[2]<0 || td.fVoxSlices[2]>fIbz-1) return 0; // outside range
624  dmin[2] = (fZb[td.fVoxSlices[2]]-point[2])*td.fVoxInvdir[2];
625  }
626  isZlimit = (dmin[2]>maxstep)?kTRUE:kFALSE;
627 // printf("---> Z : priority=%i, slice=%i/%i inc=%i\n",
628 // fPriority[2], td.fVoxSlices[2], fIbz-2, td.fVoxInc[2]);
629 // printf("4- step to next Z (%i) = %g\n", (Int_t)isZlimit, dmin[2]);
630 
631  //---> check if propagation to next slice on this axis is forced
632  if ((td.fVoxSlices[2]==-1) || (td.fVoxSlices[2]==fIbz-1)) {
633  isForcedZ = kTRUE;
634  dforced[2] = dmin[2];
635  iforced++;
636 // printf(" FORCED 1\n");
637  if (isZlimit) return 0;
638  } else {
639  if (fPriority[2]==2) {
640  // if no candidates in current slice, force next slice
641  if (fNsliceZ[td.fVoxSlices[2]]==0) {
642  isForcedZ = kTRUE;
643  dforced[2] = dmin[2];
644  iforced++;
645 // printf(" FORCED 2\n");
646  if (isZlimit) return 0;
647  }
648  }
649  }
650  } else {
651  // no slices on this axis -> bounding box limit
652 // printf(" No slice on Z\n");
653  dmin[2] = td.fVoxLimits[2];
654  isZlimit = kTRUE;
655  }
656  //---> We are done with checking. See which is the closest slice.
657  // First check if some slice is forced
658 
659  Double_t dslice = 0;
660  Int_t islice = 0;
661  if (iforced) {
662  // some slice is forced
663  if (isForcedX) {
664  // X forced
665  dslice = dforced[0];
666  islice = 0;
667  if (isForcedY) {
668  // X+Y forced
669  if (dforced[1]>dslice) {
670  dslice = dforced[1];
671  islice = 1;
672  }
673  if (isForcedZ) {
674  // X+Y+Z forced
675  if (dforced[2]>dslice) {
676  dslice = dforced[2];
677  islice = 2;
678  }
679  }
680  } else {
681  // X forced
682  if (isForcedZ) {
683  // X+Z forced
684  if (dforced[2]>dslice) {
685  dslice = dforced[2];
686  islice = 2;
687  }
688  }
689  }
690  } else {
691  if (isForcedY) {
692  // Y forced
693  dslice = dforced[1];
694  islice = 1;
695  if (isForcedZ) {
696  // Y+Z forced
697  if (dforced[2]>dslice) {
698  dslice = dforced[2];
699  islice = 2;
700  }
701  }
702  } else {
703  // Z forced
704  dslice = dforced[2];
705  islice = 2;
706  }
707  }
708  } else {
709  // Nothing forced -> get minimum distance
710  islice = TMath::LocMin(3, dmin);
711  dslice = dmin[islice];
712  if (dslice>=maxstep) {
713 // printf("DSLICE > MAXSTEP -> EXIT\n");
714  return 0;
715  }
716  }
717 // printf("5- islicenext=%i DSLICE=%g\n", islice, dslice);
718  Double_t xptnew;
719  Int_t *new_list; // list of new candidates
720  UChar_t *slice1 = 0;
721  UChar_t *slice2 = 0;
722  Int_t ndd[2] = {0,0};
723  Int_t islices = 0;
724  Bool_t left;
725  switch (islice) {
726  case 0:
727  if (isXlimit) return 0;
728  // increment/decrement X slice
729  td.fVoxSlices[0]=dind[0];
730  if (iforced) {
731  // we have to recompute Y and Z slices
732  if (dslice>td.fVoxLimits[1]) return 0;
733  if (dslice>td.fVoxLimits[2]) return 0;
734  if ((dslice>dmin[1]) && td.fVoxInc[1]) {
735  xptnew = point[1]+dslice/td.fVoxInvdir[1];
736  // printf(" recomputing Y slice, pos=%g\n", xptnew);
737  while (1) {
738  td.fVoxSlices[1] += td.fVoxInc[1];
739  if (td.fVoxInc[1]==1) {
740  if (td.fVoxSlices[1]<-1 || td.fVoxSlices[1]>fIby-2) break; // outside range
741  if (fYb[td.fVoxSlices[1]+1]>=xptnew) break;
742  } else {
743  if (td.fVoxSlices[1]<0 || td.fVoxSlices[1]>fIby-1) break; // outside range
744  if (fYb[td.fVoxSlices[1]]<= xptnew) break;
745  }
746  }
747  // printf(" %i/%i\n", td.fVoxSlices[1], fIby-2);
748  }
749  if ((dslice>dmin[2]) && td.fVoxInc[2]) {
750  xptnew = point[2]+dslice/td.fVoxInvdir[2];
751  // printf(" recomputing Z slice, pos=%g\n", xptnew);
752  while (1) {
753  td.fVoxSlices[2] += td.fVoxInc[2];
754  if (td.fVoxInc[2]==1) {
755  if (td.fVoxSlices[2]<-1 || td.fVoxSlices[2]>fIbz-2) break; // outside range
756  if (fZb[td.fVoxSlices[2]+1]>=xptnew) break;
757  } else {
758  if (td.fVoxSlices[2]<0 || td.fVoxSlices[2]>fIbz-1) break; // outside range
759  if (fZb[td.fVoxSlices[2]]<= xptnew) break;
760  }
761  }
762  // printf(" %i/%i\n", td.fVoxSlices[2], fIbz-2);
763  }
764  }
765  // new indices are set -> Get new candidates
766  if (fPriority[0]==1) {
767  // we are entering the unique slice on this axis
768  //---> intersect and store Y and Z
769  if (fPriority[1]==2) {
770  if (td.fVoxSlices[1]<0 || td.fVoxSlices[1]>=fIby-1) return td.fVoxCheckList; // outside range
771  ndd[0] = fNsliceY[td.fVoxSlices[1]];
772  if (!ndd[0]) return td.fVoxCheckList;
773  slice1 = &fIndcY[fOBy[td.fVoxSlices[1]]];
774  islices++;
775  }
776  if (fPriority[2]==2) {
777  if (td.fVoxSlices[2]<0 || td.fVoxSlices[2]>=fIbz-1) return td.fVoxCheckList; // outside range
778  ndd[1] = fNsliceZ[td.fVoxSlices[2]];
779  if (!ndd[1]) return td.fVoxCheckList;
780  islices++;
781  if (slice1) {
782  slice2 = &fIndcZ[fOBz[td.fVoxSlices[2]]];
783  } else {
784  slice1 = &fIndcZ[fOBz[td.fVoxSlices[2]]];
785  ndd[0] = ndd[1];
786  }
787  }
788  if (islices<=1) {
789  IntersectAndStore(ndd[0], slice1, td);
790  } else {
791  IntersectAndStore(ndd[0], slice1, ndd[1], slice2, td);
792  }
793  ncheck = td.fVoxNcandidates;
794  return td.fVoxCheckList;
795  }
796  // We got into a new slice -> Get only new candidates
797  left = (td.fVoxInc[0]>0)?kTRUE:kFALSE;
798  new_list = GetExtraX(td.fVoxSlices[0], left, ncheck);
799 // printf(" New list on X : %i new candidates\n", ncheck);
800  if (!ncheck) return td.fVoxCheckList;
801  if (fPriority[1]==2) {
802  if (td.fVoxSlices[1]<0 || td.fVoxSlices[1]>=fIby-1) {
803  ncheck = 0;
804  return td.fVoxCheckList; // outside range
805  }
806  ndd[0] = fNsliceY[td.fVoxSlices[1]];
807  if (!ndd[0]) {
808  ncheck = 0;
809  return td.fVoxCheckList;
810  }
811  slice1 = &fIndcY[fOBy[td.fVoxSlices[1]]];
812  islices++;
813  }
814  if (fPriority[2]==2) {
815  if (td.fVoxSlices[2]<0 || td.fVoxSlices[2]>=fIbz-1) {
816  ncheck = 0;
817  return td.fVoxCheckList; // outside range
818  }
819  ndd[1] = fNsliceZ[td.fVoxSlices[2]];
820  if (!ndd[1]) {
821  ncheck = 0;
822  return td.fVoxCheckList;
823  }
824  islices++;
825  if (slice1) {
826  slice2 = &fIndcZ[fOBz[td.fVoxSlices[2]]];
827  } else {
828  slice1 = &fIndcZ[fOBz[td.fVoxSlices[2]]];
829  ndd[0] = ndd[1];
830  }
831  }
832  if (!islices) return GetValidExtra(new_list, ncheck, td);
833  if (islices==1) {
834  return GetValidExtra(ndd[0], slice1, new_list, ncheck,td);
835  } else {
836  return GetValidExtra(ndd[0], slice1, ndd[1], slice2, new_list, ncheck, td);
837  }
838  case 1:
839  if (isYlimit) return 0;
840  // increment/decrement Y slice
841  td.fVoxSlices[1]=dind[1];
842  if (iforced) {
843  // we have to recompute X and Z slices
844  if (dslice>td.fVoxLimits[0]) return 0;
845  if (dslice>td.fVoxLimits[2]) return 0;
846  if ((dslice>dmin[0]) && td.fVoxInc[0]) {
847  xptnew = point[0]+dslice/td.fVoxInvdir[0];
848 // printf(" recomputing X slice, pos=%g\n", xptnew);
849  while (1) {
850  td.fVoxSlices[0] += td.fVoxInc[0];
851  if (td.fVoxInc[0]==1) {
852  if (td.fVoxSlices[0]<-1 || td.fVoxSlices[0]>fIbx-2) break; // outside range
853  if (fXb[td.fVoxSlices[0]+1]>=xptnew) break;
854  } else {
855  if (td.fVoxSlices[0]<0 || td.fVoxSlices[0]>fIbx-1) break; // outside range
856  if (fXb[td.fVoxSlices[0]]<= xptnew) break;
857  }
858  }
859 // printf(" %i/%i\n", td.fVoxSlices[0], fIbx-2);
860  }
861  if ((dslice>dmin[2]) && td.fVoxInc[2]) {
862  xptnew = point[2]+dslice/td.fVoxInvdir[2];
863 // printf(" recomputing Z slice, pos=%g\n", xptnew);
864  while (1) {
865  td.fVoxSlices[2] += td.fVoxInc[2];
866  if (td.fVoxInc[2]==1) {
867  if (td.fVoxSlices[2]<-1 || td.fVoxSlices[2]>fIbz-2) break; // outside range
868  if (fZb[td.fVoxSlices[2]+1]>=xptnew) break;
869  } else {
870  if (td.fVoxSlices[2]<0 || td.fVoxSlices[2]>fIbz-1) break; // outside range
871  if (fZb[td.fVoxSlices[2]]<= xptnew) break;
872  }
873  }
874 // printf(" %i/%i\n", td.fVoxSlices[2], fIbz-2);
875  }
876  }
877  // new indices are set -> Get new candidates
878  if (fPriority[1]==1) {
879  // we are entering the unique slice on this axis
880  //---> intersect and store X and Z
881  if (fPriority[0]==2) {
882  if (td.fVoxSlices[0]<0 || td.fVoxSlices[0]>=fIbx-1) return td.fVoxCheckList; // outside range
883  ndd[0] = fNsliceX[td.fVoxSlices[0]];
884  if (!ndd[0]) return td.fVoxCheckList;
885  slice1 = &fIndcX[fOBx[td.fVoxSlices[0]]];
886  islices++;
887  }
888  if (fPriority[2]==2) {
889  if (td.fVoxSlices[2]<0 || td.fVoxSlices[2]>=fIbz-1) return td.fVoxCheckList; // outside range
890  ndd[1] = fNsliceZ[td.fVoxSlices[2]];
891  if (!ndd[1]) return td.fVoxCheckList;
892  islices++;
893  if (slice1) {
894  slice2 = &fIndcZ[fOBz[td.fVoxSlices[2]]];
895  } else {
896  slice1 = &fIndcZ[fOBz[td.fVoxSlices[2]]];
897  ndd[0] = ndd[1];
898  }
899  }
900  if (islices<=1) {
901  IntersectAndStore(ndd[0], slice1, td);
902  } else {
903  IntersectAndStore(ndd[0], slice1, ndd[1], slice2, td);
904  }
905  ncheck = td.fVoxNcandidates;
906  return td.fVoxCheckList;
907  }
908  // We got into a new slice -> Get only new candidates
909  left = (td.fVoxInc[1]>0)?kTRUE:kFALSE;
910  new_list = GetExtraY(td.fVoxSlices[1], left, ncheck);
911 // printf(" New list on Y : %i new candidates\n", ncheck);
912  if (!ncheck) return td.fVoxCheckList;
913  if (fPriority[0]==2) {
914  if (td.fVoxSlices[0]<0 || td.fVoxSlices[0]>=fIbx-1) {
915  ncheck = 0;
916  return td.fVoxCheckList; // outside range
917  }
918  ndd[0] = fNsliceX[td.fVoxSlices[0]];
919  if (!ndd[0]) {
920  ncheck = 0;
921  return td.fVoxCheckList;
922  }
923  slice1 = &fIndcX[fOBx[td.fVoxSlices[0]]];
924  islices++;
925  }
926  if (fPriority[2]==2) {
927  if (td.fVoxSlices[2]<0 || td.fVoxSlices[2]>=fIbz-1) {
928  ncheck = 0;
929  return td.fVoxCheckList; // outside range
930  }
931  ndd[1] = fNsliceZ[td.fVoxSlices[2]];
932  if (!ndd[1]) {
933  ncheck = 0;
934  return td.fVoxCheckList;
935  }
936  islices++;
937  if (slice1) {
938  slice2 = &fIndcZ[fOBz[td.fVoxSlices[2]]];
939  } else {
940  slice1 = &fIndcZ[fOBz[td.fVoxSlices[2]]];
941  ndd[0] = ndd[1];
942  }
943  }
944  if (!islices) return GetValidExtra(new_list, ncheck, td);
945  if (islices==1) {
946  return GetValidExtra(ndd[0], slice1, new_list, ncheck, td);
947  } else {
948  return GetValidExtra(ndd[0], slice1, ndd[1], slice2, new_list, ncheck, td);
949  }
950  case 2:
951  if (isZlimit) return 0;
952  // increment/decrement Z slice
953  td.fVoxSlices[2]=dind[2];
954  if (iforced) {
955  // we have to recompute Y and X slices
956  if (dslice>td.fVoxLimits[1]) return 0;
957  if (dslice>td.fVoxLimits[0]) return 0;
958  if ((dslice>dmin[1]) && td.fVoxInc[1]) {
959  xptnew = point[1]+dslice/td.fVoxInvdir[1];
960  // printf(" recomputing Y slice, pos=%g\n", xptnew);
961  while (1) {
962  td.fVoxSlices[1] += td.fVoxInc[1];
963  if (td.fVoxInc[1]==1) {
964  if (td.fVoxSlices[1]<-1 || td.fVoxSlices[1]>fIby-2) break; // outside range
965  if (fYb[td.fVoxSlices[1]+1]>=xptnew) break;
966  } else {
967  if (td.fVoxSlices[1]<0 || td.fVoxSlices[1]>fIby-1) break; // outside range
968  if (fYb[td.fVoxSlices[1]]<= xptnew) break;
969  }
970  }
971  // printf(" %i/%i\n", td.fVoxSlices[1], fIby-2);
972  }
973  if ((dslice>dmin[0]) && td.fVoxInc[0]) {
974  xptnew = point[0]+dslice/td.fVoxInvdir[0];
975  // printf(" recomputing X slice, pos=%g\n", xptnew);
976  while (1) {
977  td.fVoxSlices[0] += td.fVoxInc[0];
978  if (td.fVoxInc[0]==1) {
979  if (td.fVoxSlices[0]<-1 || td.fVoxSlices[0]>fIbx-2) break; // outside range
980  if (fXb[td.fVoxSlices[0]+1]>=xptnew) break;
981  } else {
982  if (td.fVoxSlices[0]<0 || td.fVoxSlices[0]>fIbx-1) break; // outside range
983  if (fXb[td.fVoxSlices[0]]<= xptnew) break;
984  }
985  }
986 // printf(" %i/%i\n", td.fVoxSlices[0], fIbx-2);
987  }
988  }
989  // new indices are set -> Get new candidates
990  if (fPriority[2]==1) {
991  // we are entering the unique slice on this axis
992  //---> intersect and store Y and X
993  if (fPriority[1]==2) {
994  if (td.fVoxSlices[1]<0 || td.fVoxSlices[1]>=fIby-1) return td.fVoxCheckList; // outside range
995  ndd[0] = fNsliceY[td.fVoxSlices[1]];
996  if (!ndd[0]) return td.fVoxCheckList;
997  slice1 = &fIndcY[fOBy[td.fVoxSlices[1]]];
998  islices++;
999  }
1000  if (fPriority[0]==2) {
1001  if (td.fVoxSlices[0]<0 || td.fVoxSlices[0]>=fIbx-1) return td.fVoxCheckList; // outside range
1002  ndd[1] = fNsliceX[td.fVoxSlices[0]];
1003  if (!ndd[1]) return td.fVoxCheckList;
1004  islices++;
1005  if (slice1) {
1006  slice2 = &fIndcX[fOBx[td.fVoxSlices[0]]];
1007  } else {
1008  slice1 = &fIndcX[fOBx[td.fVoxSlices[0]]];
1009  ndd[0] = ndd[1];
1010  }
1011  }
1012  if (islices<=1) {
1013  IntersectAndStore(ndd[0], slice1, td);
1014  } else {
1015  IntersectAndStore(ndd[0], slice1, ndd[1], slice2, td);
1016  }
1017  ncheck = td.fVoxNcandidates;
1018  return td.fVoxCheckList;
1019  }
1020  // We got into a new slice -> Get only new candidates
1021  left = (td.fVoxInc[2]>0)?kTRUE:kFALSE;
1022  new_list = GetExtraZ(td.fVoxSlices[2], left, ncheck);
1023 // printf(" New list on Z : %i new candidates\n", ncheck);
1024  if (!ncheck) return td.fVoxCheckList;
1025  if (fPriority[1]==2) {
1026  if (td.fVoxSlices[1]<0 || td.fVoxSlices[1]>=fIby-1) {
1027  ncheck = 0;
1028  return td.fVoxCheckList; // outside range
1029  }
1030  ndd[0] = fNsliceY[td.fVoxSlices[1]];
1031  if (!ndd[0]) {
1032  ncheck = 0;
1033  return td.fVoxCheckList;
1034  }
1035  slice1 = &fIndcY[fOBy[td.fVoxSlices[1]]];
1036  islices++;
1037  }
1038  if (fPriority[0]==2) {
1039  if (td.fVoxSlices[0]<0 || td.fVoxSlices[0]>=fIbx-1) {
1040  ncheck = 0;
1041  return td.fVoxCheckList; // outside range
1042  }
1043  ndd[1] = fNsliceX[td.fVoxSlices[0]];
1044  if (!ndd[1]) {
1045  ncheck = 0;
1046  return td.fVoxCheckList;
1047  }
1048  islices++;
1049  if (slice1) {
1050  slice2 = &fIndcX[fOBx[td.fVoxSlices[0]]];
1051  } else {
1052  slice1 = &fIndcX[fOBx[td.fVoxSlices[0]]];
1053  ndd[0] = ndd[1];
1054  }
1055  }
1056  if (!islices) return GetValidExtra(new_list, ncheck, td);
1057  if (islices==1) {
1058  return GetValidExtra(ndd[0], slice1, new_list, ncheck, td);
1059  } else {
1060  return GetValidExtra(ndd[0], slice1, ndd[1], slice2, new_list, ncheck, td);
1061  }
1062  default:
1063  Error("GetNextCandidates", "Invalid islice=%i inside %s", islice, fVolume->GetName());
1064  }
1065  return 0;
1066 }
1067 
1068 ////////////////////////////////////////////////////////////////////////////////
1069 /// get the list in the next voxel crossed by a ray
1070 
1072 {
1073  if (NeedRebuild()) {
1074  TGeoVoxelFinder *vox = (TGeoVoxelFinder*)this;
1075  vox->Voxelize();
1076  fVolume->FindOverlaps();
1077  }
1078  td.fVoxCurrent = 0;
1079 // printf("###Sort crossed voxels for %s\n", fVolume->GetName());
1080  td.fVoxNcandidates = 0;
1081  Int_t loc = 1+((fVolume->GetNdaughters()-1)>>3);
1082 // printf(" LOC=%i\n", loc*sizeof(UChar_t));
1083 // UChar_t *bits = gGeoManager->GetBits();
1084  memset(td.fVoxBits1, 0, loc);
1085  memset(td.fVoxInc, 0, 3*sizeof(Int_t));
1086  for (Int_t i=0; i<3; i++) {
1087  td.fVoxInvdir[i] = TGeoShape::Big();
1088  if (TMath::Abs(dir[i])<1E-10) continue;
1089  td.fVoxInc[i] = (dir[i]>0)?1:-1;
1090  td.fVoxInvdir[i] = 1./dir[i];
1091  }
1092  Bool_t flag = GetIndices(point, td);
1093  TGeoBBox *box = (TGeoBBox*)(fVolume->GetShape());
1094  const Double_t *box_orig = box->GetOrigin();
1095  if (td.fVoxInc[0]==0) {
1096  td.fVoxLimits[0] = TGeoShape::Big();
1097  } else {
1098  if (td.fVoxSlices[0]==-2) {
1099  // no slice on this axis -> get limit to bounding box limit
1100  td.fVoxLimits[0] = (box_orig[0]-point[0]+td.fVoxInc[0]*box->GetDX())*td.fVoxInvdir[0];
1101  } else {
1102  if (td.fVoxInc[0]==1) {
1103  td.fVoxLimits[0] = (fXb[fIbx-1]-point[0])*td.fVoxInvdir[0];
1104  } else {
1105  td.fVoxLimits[0] = (fXb[0]-point[0])*td.fVoxInvdir[0];
1106  }
1107  }
1108  }
1109  if (td.fVoxInc[1]==0) {
1110  td.fVoxLimits[1] = TGeoShape::Big();
1111  } else {
1112  if (td.fVoxSlices[1]==-2) {
1113  // no slice on this axis -> get limit to bounding box limit
1114  td.fVoxLimits[1] = (box_orig[1]-point[1]+td.fVoxInc[1]*box->GetDY())*td.fVoxInvdir[1];
1115  } else {
1116  if (td.fVoxInc[1]==1) {
1117  td.fVoxLimits[1] = (fYb[fIby-1]-point[1])*td.fVoxInvdir[1];
1118  } else {
1119  td.fVoxLimits[1] = (fYb[0]-point[1])*td.fVoxInvdir[1];
1120  }
1121  }
1122  }
1123  if (td.fVoxInc[2]==0) {
1124  td.fVoxLimits[2] = TGeoShape::Big();
1125  } else {
1126  if (td.fVoxSlices[2]==-2) {
1127  // no slice on this axis -> get limit to bounding box limit
1128  td.fVoxLimits[2] = (box_orig[2]-point[2]+td.fVoxInc[2]*box->GetDZ())*td.fVoxInvdir[2];
1129  } else {
1130  if (td.fVoxInc[2]==1) {
1131  td.fVoxLimits[2] = (fZb[fIbz-1]-point[2])*td.fVoxInvdir[2];
1132  } else {
1133  td.fVoxLimits[2] = (fZb[0]-point[2])*td.fVoxInvdir[2];
1134  }
1135  }
1136  }
1137 
1138  if (!flag) {
1139 // printf(" NO candidates in first voxel\n");
1140 // printf(" bits[0]=%i\n", bits[0]);
1141  return;
1142  }
1143 // printf(" current slices : %i %i %i\n", td.fVoxSlices[0], td.fVoxSlices[1], td.fVoxSlices[2]);
1144  Int_t nd[3];
1145  Int_t islices = 0;
1146  memset(&nd[0], 0, 3*sizeof(Int_t));
1147  UChar_t *slicex = 0;
1148  if (fPriority[0]==2) {
1149  nd[0] = fNsliceX[td.fVoxSlices[0]];
1150  slicex=&fIndcX[fOBx[td.fVoxSlices[0]]];
1151  islices++;
1152  }
1153  UChar_t *slicey = 0;
1154  if (fPriority[1]==2) {
1155  nd[1] = fNsliceY[td.fVoxSlices[1]];
1156  islices++;
1157  if (slicex) {
1158  slicey=&fIndcY[fOBy[td.fVoxSlices[1]]];
1159  } else {
1160  slicex=&fIndcY[fOBy[td.fVoxSlices[1]]];
1161  nd[0] = nd[1];
1162  }
1163  }
1164  UChar_t *slicez = 0;
1165  if (fPriority[2]==2) {
1166  nd[2] = fNsliceZ[td.fVoxSlices[2]];
1167  islices++;
1168  if (slicex && slicey) {
1169  slicez=&fIndcZ[fOBz[td.fVoxSlices[2]]];
1170  } else {
1171  if (slicex) {
1172  slicey=&fIndcZ[fOBz[td.fVoxSlices[2]]];
1173  nd[1] = nd[2];
1174  } else {
1175  slicex=&fIndcZ[fOBz[td.fVoxSlices[2]]];
1176  nd[0] = nd[2];
1177  }
1178  }
1179  }
1180 // printf("Ndaughters in first voxel : %i %i %i\n", nd[0], nd[1], nd[2]);
1181  switch (islices) {
1182  case 0:
1183  Error("SortCrossedVoxels", "no slices for %s", fVolume->GetName());
1184 // printf("Slices :(%i,%i,%i) Priority:(%i,%i,%i)\n", td.fVoxSlices[0], td.fVoxSlices[1], td.fVoxSlices[2], fPriority[0], fPriority[1], fPriority[2]);
1185  return;
1186  case 1:
1187  IntersectAndStore(nd[0], slicex, td);
1188  break;
1189  case 2:
1190  IntersectAndStore(nd[0], slicex, nd[1], slicey, td);
1191  break;
1192  default:
1193  IntersectAndStore(nd[0], slicex, nd[1], slicey, nd[2], slicez, td);
1194  }
1195 // printf(" bits[0]=%i END\n", bits[0]);
1196 // if (td.fVoxNcandidates) {
1197 // printf(" candidates for first voxel :\n");
1198 // for (Int_t i=0; i<td.fVoxNcandidates; i++) printf(" %i\n", td.fVoxCheckList[i]);
1199 // }
1200 }
1201 
1202 ////////////////////////////////////////////////////////////////////////////////
1203 /// get the list of daughter indices for which point is inside their bbox
1204 
1206 {
1207  if (NeedRebuild()) {
1208  Voxelize();
1209  fVolume->FindOverlaps();
1210  }
1211  if (fVolume->GetNdaughters() == 1) {
1212  if (fXb) {
1213  if (point[0]<fXb[0] || point[0]>fXb[1]) return 0;
1214  }
1215  if (fYb) {
1216  if (point[1]<fYb[0] || point[1]>fYb[1]) return 0;
1217  }
1218 
1219  if (fZb) {
1220  if (point[2]<fZb[0] || point[2]>fZb[1]) return 0;
1221  }
1222  td.fVoxCheckList[0] = 0;
1223  nelem = 1;
1224  return td.fVoxCheckList;
1225  }
1226  Int_t nslices = 0;
1227  UChar_t *slice1 = 0;
1228  UChar_t *slice2 = 0;
1229  UChar_t *slice3 = 0;
1230  Int_t nd[3] = {0,0,0};
1231  Int_t im;
1232  if (fPriority[0]) {
1233  im = TMath::BinarySearch(fIbx, fXb, point[0]);
1234  if ((im==-1) || (im==fIbx-1)) return 0;
1235  if (fPriority[0]==2) {
1236  nd[0] = fNsliceX[im];
1237  if (!nd[0]) return 0;
1238  nslices++;
1239  slice1 = &fIndcX[fOBx[im]];
1240  }
1241  }
1242 
1243  if (fPriority[1]) {
1244  im = TMath::BinarySearch(fIby, fYb, point[1]);
1245  if ((im==-1) || (im==fIby-1)) return 0;
1246  if (fPriority[1]==2) {
1247  nd[1] = fNsliceY[im];
1248  if (!nd[1]) return 0;
1249  nslices++;
1250  if (slice1) {
1251  slice2 = &fIndcY[fOBy[im]];
1252  } else {
1253  slice1 = &fIndcY[fOBy[im]];
1254  nd[0] = nd[1];
1255  }
1256  }
1257  }
1258 
1259  if (fPriority[2]) {
1260  im = TMath::BinarySearch(fIbz, fZb, point[2]);
1261  if ((im==-1) || (im==fIbz-1)) return 0;
1262  if (fPriority[2]==2) {
1263  nd[2] = fNsliceZ[im];
1264  if (!nd[2]) return 0;
1265  nslices++;
1266  if (slice1 && slice2) {
1267  slice3 = &fIndcZ[fOBz[im]];
1268  } else {
1269  if (slice1) {
1270  slice2 = &fIndcZ[fOBz[im]];
1271  nd[1] = nd[2];
1272  } else {
1273  slice1 = &fIndcZ[fOBz[im]];
1274  nd[0] = nd[2];
1275  }
1276  }
1277  }
1278  }
1279  nelem = 0;
1280 // Int_t i = 0;
1281  Bool_t intersect = kFALSE;
1282  switch (nslices) {
1283  case 0:
1284  Error("GetCheckList", "No slices for %s", fVolume->GetName());
1285  return 0;
1286  case 1:
1287  intersect = Intersect(nd[0], slice1, nelem, td.fVoxCheckList);
1288  break;
1289  case 2:
1290  intersect = Intersect(nd[0], slice1, nd[1], slice2, nelem, td.fVoxCheckList);
1291  break;
1292  default:
1293  intersect = Intersect(nd[0], slice1, nd[1], slice2, nd[2], slice3, nelem, td.fVoxCheckList);
1294  }
1295  if (intersect) return td.fVoxCheckList;
1296  return 0;
1297 }
1298 
1299 ////////////////////////////////////////////////////////////////////////////////
1300 /// get the list of candidates in voxel (i,j,k) - no check
1301 
1303 {
1304  UChar_t *slice1 = 0;
1305  UChar_t *slice2 = 0;
1306  UChar_t *slice3 = 0;
1307  Int_t nd[3] = {0,0,0};
1308  Int_t nslices = 0;
1309  if (fPriority[0]==2) {
1310  nd[0] = fNsliceX[i];
1311  if (!nd[0]) return 0;
1312  nslices++;
1313  slice1 = &fIndcX[fOBx[i]];
1314  }
1315 
1316  if (fPriority[1]==2) {
1317  nd[1] = fNsliceY[j];
1318  if (!nd[1]) return 0;
1319  nslices++;
1320  if (slice1) {
1321  slice2 = &fIndcY[fOBy[j]];
1322  } else {
1323  slice1 = &fIndcY[fOBy[j]];
1324  nd[0] = nd[1];
1325  }
1326  }
1327 
1328  if (fPriority[2]==2) {
1329  nd[2] = fNsliceZ[k];
1330  if (!nd[2]) return 0;
1331  nslices++;
1332  if (slice1 && slice2) {
1333  slice3 = &fIndcZ[fOBz[k]];
1334  } else {
1335  if (slice1) {
1336  slice2 = &fIndcZ[fOBz[k]];
1337  nd[1] = nd[2];
1338  } else {
1339  slice1 = &fIndcZ[fOBz[k]];
1340  nd[0] = nd[2];
1341  }
1342  }
1343  }
1344  Bool_t intersect = kFALSE;
1345  switch (nslices) {
1346  case 0:
1347  Error("GetCheckList", "No slices for %s", fVolume->GetName());
1348  return 0;
1349  case 1:
1350  intersect = Intersect(nd[0], slice1, ncheck, td.fVoxCheckList);
1351  break;
1352  case 2:
1353  intersect = Intersect(nd[0], slice1, nd[1], slice2, ncheck, td.fVoxCheckList);
1354  break;
1355  default:
1356  intersect = Intersect(nd[0], slice1, nd[1], slice2, nd[2], slice3, ncheck, td.fVoxCheckList);
1357  }
1358  if (intersect) return td.fVoxCheckList;
1359  return 0;
1360 }
1361 
1362 ////////////////////////////////////////////////////////////////////////////////
1363 /// get the list of new candidates for the next voxel crossed by current ray
1364 /// printf("### GetNextVoxel\n");
1365 
1366 Int_t *TGeoVoxelFinder::GetNextVoxel(const Double_t *point, const Double_t * /*dir*/, Int_t &ncheck, TGeoStateInfo &td)
1367 {
1368  if (NeedRebuild()) {
1369  Voxelize();
1370  fVolume->FindOverlaps();
1371  }
1372  if (td.fVoxCurrent==0) {
1373 // printf(">>> first voxel, %i candidates\n", ncheck);
1374 // printf(" bits[0]=%i\n", gGeoManager->GetBits()[0]);
1375  td.fVoxCurrent++;
1376  ncheck = td.fVoxNcandidates;
1377  return td.fVoxCheckList;
1378  }
1379  td.fVoxCurrent++;
1380 // printf(">>> voxel %i\n", td.fCurrentVoxel);
1381  // Get slices for next voxel
1382 // printf("before - td.fSlices : %i %i %i\n", td.fSlices[0], td.fSlices[1], td.fSlices[2]);
1383  return GetNextCandidates(point, ncheck, td);
1384 }
1385 
1386 ////////////////////////////////////////////////////////////////////////////////
1387 /// return the list of nodes corresponding to one array of bits
1388 
1390 {
1391  Int_t nd = fVolume->GetNdaughters(); // also number of bits to scan
1392  nf = 0;
1393  Int_t nbytes = 1+((nd-1)>>3);
1394  Int_t current_byte;
1395  Int_t current_bit;
1396  UChar_t byte;
1397  Bool_t ibreak = kFALSE;
1398  for (current_byte=0; current_byte<nbytes; current_byte++) {
1399  byte = array1[current_byte];
1400  if (!byte) continue;
1401  for (current_bit=0; current_bit<8; current_bit++) {
1402  if (byte & (1<<current_bit)) {
1403  result[nf++] = (current_byte<<3)+current_bit;
1404  if (nf==n1) {
1405  ibreak = kTRUE;
1406  break;
1407  }
1408  }
1409  }
1410  if (ibreak) return kTRUE;
1411  }
1412  return kTRUE;
1413 }
1414 
1415 ////////////////////////////////////////////////////////////////////////////////
1416 /// return the list of nodes corresponding to one array of bits
1417 
1419 {
1420  Int_t nd = fVolume->GetNdaughters(); // also number of bits to scan
1421 // UChar_t *bits = gGeoManager->GetBits();
1422  td.fVoxNcandidates = 0;
1423  Int_t nbytes = 1+((nd-1)>>3);
1424  if (!array1) {
1425  memset(td.fVoxBits1, 0xFF, nbytes*sizeof(UChar_t));
1426  while (td.fVoxNcandidates<nd) {
1428  ++td.fVoxNcandidates;
1429  }
1430  return kTRUE;
1431  }
1432  memcpy(td.fVoxBits1, array1, nbytes*sizeof(UChar_t));
1433  Int_t current_byte;
1434  Int_t current_bit;
1435  UChar_t byte;
1436  Bool_t ibreak = kFALSE;
1437  Int_t icand;
1438  for (current_byte=0; current_byte<nbytes; current_byte++) {
1439  byte = array1[current_byte];
1440  icand = current_byte<<3;
1441  if (!byte) continue;
1442  for (current_bit=0; current_bit<8; current_bit++) {
1443  if (byte & (1<<current_bit)) {
1444  td.fVoxCheckList[td.fVoxNcandidates++] = icand+current_bit;
1445  if (td.fVoxNcandidates==n1) {
1446  ibreak = kTRUE;
1447  break;
1448  }
1449  }
1450  }
1451  if (ibreak) return kTRUE;
1452  }
1453  return kTRUE;
1454 }
1455 
1456 ////////////////////////////////////////////////////////////////////////////////
1457 /// make union of older bits with new array
1458 /// printf("Union - one slice\n");
1459 
1461 {
1462  Int_t nd = fVolume->GetNdaughters(); // also number of bits to scan
1463 // UChar_t *bits = gGeoManager->GetBits();
1464  td.fVoxNcandidates = 0;
1465  Int_t nbytes = 1+((nd-1)>>3);
1466  Int_t current_byte;
1467  Int_t current_bit;
1468  UChar_t byte;
1469  Bool_t ibreak = kFALSE;
1470  for (current_byte=0; current_byte<nbytes; current_byte++) {
1471 // printf(" byte %i : bits=%i array=%i\n", current_byte, bits[current_byte], array1[current_byte]);
1472  byte = (~td.fVoxBits1[current_byte]) & array1[current_byte];
1473  if (!byte) continue;
1474  for (current_bit=0; current_bit<8; current_bit++) {
1475  if (byte & (1<<current_bit)) {
1476  td.fVoxCheckList[td.fVoxNcandidates++] = (current_byte<<3)+current_bit;
1477  if (td.fVoxNcandidates==n1) {
1478  ibreak = kTRUE;
1479  break;
1480  }
1481  }
1482  }
1483  td.fVoxBits1[current_byte] |= byte;
1484  if (ibreak) return kTRUE;
1485  }
1486  return (td.fVoxNcandidates>0);
1487 }
1488 
1489 ////////////////////////////////////////////////////////////////////////////////
1490 /// make union of older bits with new array
1491 /// printf("Union - two slices\n");
1492 
1493 Bool_t TGeoVoxelFinder::Union(Int_t /*n1*/, UChar_t *array1, Int_t /*n2*/, UChar_t *array2, TGeoStateInfo &td)
1494 {
1495  Int_t nd = fVolume->GetNdaughters(); // also number of bits to scan
1496 // UChar_t *bits = gGeoManager->GetBits();
1497  td.fVoxNcandidates = 0;
1498  Int_t nbytes = 1+((nd-1)>>3);
1499  Int_t current_byte;
1500  Int_t current_bit;
1501  UChar_t byte;
1502  for (current_byte=0; current_byte<nbytes; current_byte++) {
1503  byte = (~td.fVoxBits1[current_byte]) & (array1[current_byte] & array2[current_byte]);
1504  if (!byte) continue;
1505  for (current_bit=0; current_bit<8; current_bit++) {
1506  if (byte & (1<<current_bit)) {
1507  td.fVoxCheckList[td.fVoxNcandidates++] = (current_byte<<3)+current_bit;
1508  }
1509  }
1510  td.fVoxBits1[current_byte] |= byte;
1511  }
1512  return (td.fVoxNcandidates>0);
1513 }
1514 
1515 ////////////////////////////////////////////////////////////////////////////////
1516 /// make union of older bits with new array
1517 /// printf("Union - three slices\n");
1518 /// printf("n1=%i n2=%i n3=%i\n", n1,n2,n3);
1519 
1520 Bool_t TGeoVoxelFinder::Union(Int_t /*n1*/, UChar_t *array1, Int_t /*n2*/, UChar_t *array2, Int_t /*n3*/, UChar_t *array3, TGeoStateInfo &td)
1521 {
1522  Int_t nd = fVolume->GetNdaughters(); // also number of bits to scan
1523 // UChar_t *bits = gGeoManager->GetBits();
1524  td.fVoxNcandidates = 0;
1525  Int_t nbytes = 1+((nd-1)>>3);
1526  Int_t current_byte;
1527  Int_t current_bit;
1528  UChar_t byte;
1529  for (current_byte=0; current_byte<nbytes; current_byte++) {
1530  byte = (~td.fVoxBits1[current_byte]) & (array1[current_byte] & array2[current_byte] & array3[current_byte]);
1531  if (!byte) continue;
1532  for (current_bit=0; current_bit<8; current_bit++) {
1533  if (byte & (1<<current_bit)) {
1534  td.fVoxCheckList[td.fVoxNcandidates++] = (current_byte<<3)+current_bit;
1535  }
1536  }
1537  td.fVoxBits1[current_byte] |= byte;
1538  }
1539  return (td.fVoxNcandidates>0);
1540 }
1541 
1542 ////////////////////////////////////////////////////////////////////////////////
1543 /// return the list of nodes corresponding to the intersection of two arrays of bits
1544 
1545 Bool_t TGeoVoxelFinder::Intersect(Int_t n1, UChar_t *array1, Int_t n2, UChar_t *array2, Int_t &nf, Int_t *result)
1546 {
1547  Int_t nd = fVolume->GetNdaughters(); // also number of bits to scan
1548  nf = 0;
1549  Int_t nbytes = 1+((nd-1)>>3);
1550  Int_t current_byte;
1551  Int_t current_bit;
1552  UChar_t byte;
1553  Bool_t ibreak = kFALSE;
1554  for (current_byte=0; current_byte<nbytes; current_byte++) {
1555  byte = array1[current_byte] & array2[current_byte];
1556  if (!byte) continue;
1557  for (current_bit=0; current_bit<8; current_bit++) {
1558  if (byte & (1<<current_bit)) {
1559  result[nf++] = (current_byte<<3)+current_bit;
1560  if ((nf==n1) || (nf==n2)) {
1561  ibreak = kTRUE;
1562  break;
1563  }
1564  }
1565  }
1566  if (ibreak) return kTRUE;
1567  }
1568  return (nf>0);
1569 }
1570 
1571 ////////////////////////////////////////////////////////////////////////////////
1572 /// return the list of nodes corresponding to the intersection of two arrays of bits
1573 
1575 {
1576  Int_t nd = fVolume->GetNdaughters(); // also number of bits to scan
1577 // UChar_t *bits = gGeoManager->GetBits();
1578  td.fVoxNcandidates = 0;
1579  Int_t nbytes = 1+((nd-1)>>3);
1580 // memset(bits, 0, nbytes*sizeof(UChar_t));
1581  Int_t current_byte;
1582  Int_t current_bit;
1583  Int_t icand;
1584  UChar_t byte;
1585  for (current_byte=0; current_byte<nbytes; current_byte++) {
1586  byte = array1[current_byte] & array2[current_byte];
1587  icand = current_byte<<3;
1588  td.fVoxBits1[current_byte] = byte;
1589  if (!byte) continue;
1590  for (current_bit=0; current_bit<8; current_bit++) {
1591  if (byte & (1<<current_bit)) {
1592  td.fVoxCheckList[td.fVoxNcandidates++] = icand+current_bit;
1593  }
1594  }
1595  }
1596  return (td.fVoxNcandidates>0);
1597 }
1598 
1599 ////////////////////////////////////////////////////////////////////////////////
1600 /// return the list of nodes corresponding to the intersection of three arrays of bits
1601 
1602 Bool_t TGeoVoxelFinder::Intersect(Int_t n1, UChar_t *array1, Int_t n2, UChar_t *array2, Int_t n3, UChar_t *array3, Int_t &nf, Int_t *result)
1603 {
1604  Int_t nd = fVolume->GetNdaughters(); // also number of bits to scan
1605  nf = 0;
1606  Int_t nbytes = 1+((nd-1)>>3);
1607  Int_t current_byte;
1608  Int_t current_bit;
1609  UChar_t byte;
1610  Bool_t ibreak = kFALSE;
1611  for (current_byte=0; current_byte<nbytes; current_byte++) {
1612  byte = array1[current_byte] & array2[current_byte] & array3[current_byte];
1613  if (!byte) continue;
1614  for (current_bit=0; current_bit<8; current_bit++) {
1615  if (byte & (1<<current_bit)) {
1616  result[nf++] = (current_byte<<3)+current_bit;
1617  if ((nf==n1) || (nf==n2) || (nf==n3)) {
1618  ibreak = kTRUE;
1619  break;
1620  }
1621  }
1622  }
1623  if (ibreak) return kTRUE;
1624  }
1625  return (nf>0);
1626 }
1627 
1628 ////////////////////////////////////////////////////////////////////////////////
1629 /// return the list of nodes corresponding to the intersection of three arrays of bits
1630 
1631 Bool_t TGeoVoxelFinder::IntersectAndStore(Int_t /*n1*/, UChar_t *array1, Int_t /*n2*/, UChar_t *array2, Int_t /*n3*/, UChar_t *array3, TGeoStateInfo &td)
1632 {
1633  Int_t nd = fVolume->GetNdaughters(); // also number of bits to scan
1634 // UChar_t *bits = gGeoManager->GetBits();
1635  td.fVoxNcandidates = 0;
1636  Int_t nbytes = 1+((nd-1)>>3);
1637 // memset(bits, 0, nbytes*sizeof(UChar_t));
1638  Int_t current_byte;
1639  Int_t current_bit;
1640  Int_t icand;
1641  UChar_t byte;
1642  for (current_byte=0; current_byte<nbytes; current_byte++) {
1643  byte = array1[current_byte] & array2[current_byte] & array3[current_byte];
1644  icand = current_byte<<3;
1645  td.fVoxBits1[current_byte] = byte;
1646  if (!byte) continue;
1647  for (current_bit=0; current_bit<8; current_bit++) {
1648  if (byte & (1<<current_bit)) {
1649  td.fVoxCheckList[td.fVoxNcandidates++] = icand+current_bit;
1650  }
1651  }
1652  }
1653  return (td.fVoxNcandidates>0);
1654 }
1655 ////////////////////////////////////////////////////////////////////////////////
1656 /// order bounding boxes along x, y, z
1657 
1659 {
1660  Int_t nd = fVolume->GetNdaughters();
1661  Int_t nperslice = 1+(nd-1)/(8*sizeof(UChar_t)); /*Nbytes per slice*/
1662  Int_t nmaxslices = 2*nd+1; // max number of slices on each axis
1663  Double_t xmin, xmax, ymin, ymax, zmin, zmax;
1664  TGeoBBox *box = (TGeoBBox*)fVolume->GetShape(); // bounding box for volume
1665  // compute range on X, Y, Z according to volume bounding box
1666  xmin = (box->GetOrigin())[0] - box->GetDX();
1667  xmax = (box->GetOrigin())[0] + box->GetDX();
1668  ymin = (box->GetOrigin())[1] - box->GetDY();
1669  ymax = (box->GetOrigin())[1] + box->GetDY();
1670  zmin = (box->GetOrigin())[2] - box->GetDZ();
1671  zmax = (box->GetOrigin())[2] + box->GetDZ();
1672  if ((xmin>=xmax) || (ymin>=ymax) || (zmin>=zmax)) {
1673  Error("SortAll", "Wrong bounding box for volume %s", fVolume->GetName());
1674  return;
1675  }
1676  Int_t id;
1677  // compute boundaries coordinates on X,Y,Z
1678  Double_t *boundaries = new Double_t[6*nd]; // list of different boundaries
1679  for (id=0; id<nd; id++) {
1680  // x boundaries
1681  boundaries[2*id] = fBoxes[6*id+3]-fBoxes[6*id];
1682  boundaries[2*id+1] = fBoxes[6*id+3]+fBoxes[6*id];
1683  // y boundaries
1684  boundaries[2*id+2*nd] = fBoxes[6*id+4]-fBoxes[6*id+1];
1685  boundaries[2*id+2*nd+1] = fBoxes[6*id+4]+fBoxes[6*id+1];
1686  // z boundaries
1687  boundaries[2*id+4*nd] = fBoxes[6*id+5]-fBoxes[6*id+2];
1688  boundaries[2*id+4*nd+1] = fBoxes[6*id+5]+fBoxes[6*id+2];
1689  }
1690  Int_t *index = new Int_t[2*nd]; // indexes for sorted boundaries on one axis
1691  UChar_t *ind = new UChar_t[nmaxslices*nperslice]; // ind[fOBx[i]] = ndghts in slice fInd[i]--fInd[i+1]
1692  Int_t *extra = new Int_t[nmaxslices*4];
1693  // extra[fOEx[i]] = nextra_to_left (i/i-1)
1694  // extra[fOEx[i]+1] = nextra_to_right (i/i+1)
1695  // Int_t *extra_to_left = extra[fOEx[i]+2]
1696  // Int_t *extra_to_right = extra[fOEx[i]+2+nextra_to_left]
1697  Double_t *temp = new Double_t[2*nd]; // temporary array to store sorted boundary positions
1698  Int_t current = 0;
1699  Int_t indextra = 0;
1700  Int_t nleft, nright;
1701  Int_t *extra_left = new Int_t[nd];
1702  Int_t *extra_right = new Int_t[nd];
1703  Double_t xxmin, xxmax, xbmin, xbmax, ddx1, ddx2;
1704  UChar_t *bits;
1705  UInt_t loc, bitnumber;
1706  UChar_t bit;
1707 
1708  // sort x boundaries
1709  Int_t ib = 0; // total number of DIFFERENT boundaries
1710  TMath::Sort(2*nd, &boundaries[0], &index[0], kFALSE);
1711  // compact common boundaries
1712  for (id=0; id<2*nd; id++) {
1713  if (!ib) {temp[ib++] = boundaries[index[id]]; continue;}
1714  if (TMath::Abs(temp[ib-1]-boundaries[index[id]])>1E-10)
1715  temp[ib++] = boundaries[index[id]];
1716  }
1717  // now find priority
1718  if (ib < 2) {
1719  Error("SortAll", "Cannot voxelize %s :less than 2 boundaries on X", fVolume->GetName());
1720  delete [] boundaries;
1721  delete [] index;
1722  delete [] ind;
1723  delete [] temp;
1724  delete [] extra;
1725  delete [] extra_left;
1726  delete [] extra_right;
1727  SetInvalid();
1728  return;
1729  }
1730  if (ib == 2) {
1731  // check range
1732  if (((temp[0]-xmin)<1E-10) && ((temp[1]-xmax)>-1E-10)) {
1733  // ordering on this axis makes no sense. Clear all arrays.
1734  fPriority[0] = 0; // always skip this axis
1735  if (fIndcX) delete [] fIndcX;
1736  fIndcX = 0;
1737  fNx = 0;
1738  if (fXb) delete [] fXb;
1739  fXb = 0;
1740  fIbx = 0;
1741  if (fOBx) delete [] fOBx;
1742  fOBx = 0;
1743  fNox = 0;
1744  } else {
1745  fPriority[0] = 1; // all in one slice
1746  }
1747  } else {
1748  fPriority[0] = 2; // check all
1749  }
1750  // store compacted boundaries
1751  if (fPriority[0]) {
1752  if (fXb) delete [] fXb;
1753  fXb = new Double_t[ib];
1754  memcpy(fXb, &temp[0], ib*sizeof(Double_t));
1755  fIbx = ib;
1756  }
1757 
1758  //--- now build the lists of nodes in each slice of X axis
1759  if (fPriority[0]==2) {
1760  memset(ind, 0, (nmaxslices*nperslice)*sizeof(UChar_t));
1761  if (fOBx) delete [] fOBx;
1762  fNox = fIbx-1; // number of different slices
1763  fOBx = new Int_t[fNox]; // offsets in ind
1764  if (fOEx) delete [] fOEx;
1765  fOEx = new Int_t[fNox]; // offsets in extra
1766  if (fNsliceX) delete [] fNsliceX;
1767  fNsliceX = new Int_t[fNox];
1768  current = 0;
1769  indextra = 0;
1770  //--- now loop all slices
1771  for (id=0; id<fNox; id++) {
1772  fOBx[id] = current; // offset in dght list for this slice
1773  fOEx[id] = indextra; // offset in exta list for this slice
1774  fNsliceX[id] = 0; // ndght in this slice
1775  extra[indextra] = extra[indextra+1] = 0; // no extra left/right
1776  nleft = nright = 0;
1777  bits = &ind[current]; // adress of bits for this slice
1778  xxmin = fXb[id];
1779  xxmax = fXb[id+1];
1780  for (Int_t ic=0; ic<nd; ic++) {
1781  xbmin = fBoxes[6*ic+3]-fBoxes[6*ic];
1782  xbmax = fBoxes[6*ic+3]+fBoxes[6*ic];
1783  ddx1 = xbmin-xxmax;
1784  if (ddx1>-1E-10) continue;
1785  ddx2 = xbmax-xxmin;
1786  if (ddx2<1E-10) continue;
1787  // daughter ic in interval
1788  //---> set the bit
1789  fNsliceX[id]++;
1790  bitnumber = (UInt_t)ic;
1791  loc = bitnumber/8;
1792  bit = bitnumber%8;
1793  bits[loc] |= 1<<bit;
1794  //---> chech if it is extra to left/right
1795  //--- left
1796  ddx1 = xbmin-xxmin;
1797  ddx2 = xbmax-xxmax;
1798  if ((id==0) || (ddx1>-1E-10)) {
1799  extra_left[nleft++] = ic;
1800  }
1801  //---right
1802  if ((id==(fNoz-1)) || (ddx2<1E-10)) {
1803  extra_right[nright++] = ic;
1804  }
1805  }
1806  //--- compute offset of next slice
1807  if (fNsliceX[id]>0) current += nperslice;
1808  //--- copy extra candidates
1809  extra[indextra] = nleft;
1810  extra[indextra+1] = nright;
1811  if (nleft) memcpy(&extra[indextra+2], extra_left, nleft*sizeof(Int_t));
1812  if (nright) memcpy(&extra[indextra+2+nleft], extra_right, nright*sizeof(Int_t));
1813  indextra += 2+nleft+nright;
1814  }
1815  if (fIndcX) delete [] fIndcX;
1816  fNx = current;
1817  fIndcX = new UChar_t[current];
1818  memcpy(fIndcX, ind, current*sizeof(UChar_t));
1819  if (fExtraX) delete [] fExtraX;
1820  fNex = indextra;
1821  if (indextra>nmaxslices*4) printf("Woops!!!\n");
1822  fExtraX = new Int_t[indextra];
1823  memcpy(fExtraX, extra, indextra*sizeof(Int_t));
1824  }
1825 
1826  // sort y boundaries
1827  ib = 0;
1828  TMath::Sort(2*nd, &boundaries[2*nd], &index[0], kFALSE);
1829  // compact common boundaries
1830  for (id=0; id<2*nd; id++) {
1831  if (!ib) {temp[ib++] = boundaries[2*nd+index[id]]; continue;}
1832  if (TMath::Abs(temp[ib-1]-boundaries[2*nd+index[id]])>1E-10)
1833  temp[ib++]=boundaries[2*nd+index[id]];
1834  }
1835  // now find priority on Y
1836  if (ib < 2) {
1837  Error("SortAll", "Cannot voxelize %s :less than 2 boundaries on Y", fVolume->GetName());
1838  delete [] boundaries;
1839  delete [] index;
1840  delete [] ind;
1841  delete [] temp;
1842  delete [] extra;
1843  delete [] extra_left;
1844  delete [] extra_right;
1845  SetInvalid();
1846  return;
1847  }
1848  if (ib == 2) {
1849  // check range
1850  if (((temp[0]-ymin)<1E-10) && ((temp[1]-ymax)>-1E-10)) {
1851  // ordering on this axis makes no sense. Clear all arrays.
1852  fPriority[1] = 0; // always skip this axis
1853  if (fIndcY) delete [] fIndcY;
1854  fIndcY = 0;
1855  fNy = 0;
1856  if (fYb) delete [] fYb;
1857  fYb = 0;
1858  fIby = 0;
1859  if (fOBy) delete [] fOBy;
1860  fOBy = 0;
1861  fNoy = 0;
1862  } else {
1863  fPriority[1] = 1; // all in one slice
1864  }
1865  } else {
1866  fPriority[1] = 2; // check all
1867  }
1868 
1869  // store compacted boundaries
1870  if (fPriority[1]) {
1871  if (fYb) delete [] fYb;
1872  fYb = new Double_t[ib];
1873  memcpy(fYb, &temp[0], ib*sizeof(Double_t));
1874  fIby = ib;
1875  }
1876 
1877 
1878  if (fPriority[1]==2) {
1879  //--- now build the lists of nodes in each slice of Y axis
1880  memset(ind, 0, (nmaxslices*nperslice)*sizeof(UChar_t));
1881  if (fOBy) delete [] fOBy;
1882  fNoy = fIby-1; // number of different slices
1883  fOBy = new Int_t[fNoy]; // offsets in ind
1884  if (fOEy) delete [] fOEy;
1885  fOEy = new Int_t[fNoy]; // offsets in extra
1886  if (fNsliceY) delete [] fNsliceY;
1887  fNsliceY = new Int_t[fNoy];
1888  current = 0;
1889  indextra = 0;
1890  //--- now loop all slices
1891  for (id=0; id<fNoy; id++) {
1892  fOBy[id] = current; // offset of dght list
1893  fOEy[id] = indextra; // offset in exta list for this slice
1894  fNsliceY[id] = 0; // ndght in this slice
1895  extra[indextra] = extra[indextra+1] = 0; // no extra left/right
1896  nleft = nright = 0;
1897  bits = &ind[current]; // adress of bits for this slice
1898  xxmin = fYb[id];
1899  xxmax = fYb[id+1];
1900  for (Int_t ic=0; ic<nd; ic++) {
1901  xbmin = fBoxes[6*ic+4]-fBoxes[6*ic+1];
1902  xbmax = fBoxes[6*ic+4]+fBoxes[6*ic+1];
1903  ddx1 = xbmin-xxmax;
1904  if (ddx1>-1E-10) continue;
1905  ddx2 = xbmax-xxmin;
1906  if (ddx2<1E-10) continue;
1907  // daughter ic in interval
1908  //---> set the bit
1909  fNsliceY[id]++;
1910  bitnumber = (UInt_t)ic;
1911  loc = bitnumber/8;
1912  bit = bitnumber%8;
1913  bits[loc] |= 1<<bit;
1914  //---> check if it is extra to left/right
1915  //--- left
1916  ddx1 = xbmin-xxmin;
1917  ddx2 = xbmax-xxmax;
1918  if ((id==0) || (ddx1>-1E-10)) {
1919  extra_left[nleft++] = ic;
1920  }
1921  //---right
1922  if ((id==(fNoz-1)) || (ddx2<1E-10)) {
1923  extra_right[nright++] = ic;
1924  }
1925  }
1926  //--- compute offset of next slice
1927  if (fNsliceY[id]>0) current += nperslice;
1928  //--- copy extra candidates
1929  extra[indextra] = nleft;
1930  extra[indextra+1] = nright;
1931  if (nleft) memcpy(&extra[indextra+2], extra_left, nleft*sizeof(Int_t));
1932  if (nright) memcpy(&extra[indextra+2+nleft], extra_right, nright*sizeof(Int_t));
1933  indextra += 2+nleft+nright;
1934  }
1935  if (fIndcY) delete [] fIndcY;
1936  fNy = current;
1937  fIndcY = new UChar_t[current];
1938  memcpy(fIndcY, &ind[0], current*sizeof(UChar_t));
1939  if (fExtraY) delete [] fExtraY;
1940  fNey = indextra;
1941  if (indextra>nmaxslices*4) printf("Woops!!!\n");
1942  fExtraY = new Int_t[indextra];
1943  memcpy(fExtraY, extra, indextra*sizeof(Int_t));
1944  }
1945 
1946  // sort z boundaries
1947  ib = 0;
1948  TMath::Sort(2*nd, &boundaries[4*nd], &index[0], kFALSE);
1949  // compact common boundaries
1950  for (id=0; id<2*nd; id++) {
1951  if (!ib) {temp[ib++] = boundaries[4*nd+index[id]]; continue;}
1952  if ((TMath::Abs(temp[ib-1]-boundaries[4*nd+index[id]]))>1E-10)
1953  temp[ib++]=boundaries[4*nd+index[id]];
1954  }
1955  // now find priority on Z
1956  if (ib < 2) {
1957  Error("SortAll", "Cannot voxelize %s :less than 2 boundaries on Z", fVolume->GetName());
1958  delete [] boundaries;
1959  delete [] index;
1960  delete [] ind;
1961  delete [] temp;
1962  delete [] extra;
1963  delete [] extra_left;
1964  delete [] extra_right;
1965  SetInvalid();
1966  return;
1967  }
1968  if (ib == 2) {
1969  // check range
1970  if (((temp[0]-zmin)<1E-10) && ((temp[1]-zmax)>-1E-10)) {
1971  // ordering on this axis makes no sense. Clear all arrays.
1972  fPriority[2] = 0;
1973  if (fIndcZ) delete [] fIndcZ;
1974  fIndcZ = 0;
1975  fNz = 0;
1976  if (fZb) delete [] fZb;
1977  fZb = 0;
1978  fIbz = 0;
1979  if (fOBz) delete [] fOBz;
1980  fOBz = 0;
1981  fNoz = 0;
1982  } else {
1983  fPriority[2] = 1; // all in one slice
1984  }
1985  } else {
1986  fPriority[2] = 2; // check all
1987  }
1988 
1989  // store compacted boundaries
1990  if (fPriority[2]) {
1991  if (fZb) delete [] fZb;
1992  fZb = new Double_t[ib];
1993  memcpy(fZb, &temp[0], ib*sizeof(Double_t));
1994  fIbz = ib;
1995  }
1996 
1997 
1998  if (fPriority[2]==2) {
1999  //--- now build the lists of nodes in each slice of Y axis
2000  memset(ind, 0, (nmaxslices*nperslice)*sizeof(UChar_t));
2001  if (fOBz) delete [] fOBz;
2002  fNoz = fIbz-1; // number of different slices
2003  fOBz = new Int_t[fNoz]; // offsets in ind
2004  if (fOEz) delete [] fOEz;
2005  fOEz = new Int_t[fNoz]; // offsets in extra
2006  if (fNsliceZ) delete [] fNsliceZ;
2007  fNsliceZ = new Int_t[fNoz];
2008  current = 0;
2009  indextra = 0;
2010  //--- now loop all slices
2011  for (id=0; id<fNoz; id++) {
2012  fOBz[id] = current; // offset of dght list
2013  fOEz[id] = indextra; // offset in exta list for this slice
2014  fNsliceZ[id] = 0; // ndght in this slice
2015  extra[indextra] = extra[indextra+1] = 0; // no extra left/right
2016  nleft = nright = 0;
2017  bits = &ind[current]; // adress of bits for this slice
2018  xxmin = fZb[id];
2019  xxmax = fZb[id+1];
2020  for (Int_t ic=0; ic<nd; ic++) {
2021  xbmin = fBoxes[6*ic+5]-fBoxes[6*ic+2];
2022  xbmax = fBoxes[6*ic+5]+fBoxes[6*ic+2];
2023  ddx1 = xbmin-xxmax;
2024  if (ddx1>-1E-10) continue;
2025  ddx2 = xbmax-xxmin;
2026  if (ddx2<1E-10) continue;
2027  // daughter ic in interval
2028  //---> set the bit
2029  fNsliceZ[id]++;
2030  bitnumber = (UInt_t)ic;
2031  loc = bitnumber/8;
2032  bit = bitnumber%8;
2033  bits[loc] |= 1<<bit;
2034  //---> check if it is extra to left/right
2035  //--- left
2036  ddx1 = xbmin-xxmin;
2037  ddx2 = xbmax-xxmax;
2038  if ((id==0) || (ddx1>-1E-10)) {
2039  extra_left[nleft++] = ic;
2040  }
2041  //---right
2042  if ((id==(fNoz-1)) || (ddx2<1E-10)) {
2043  extra_right[nright++] = ic;
2044  }
2045  }
2046  //--- compute offset of next slice
2047  if (fNsliceZ[id]>0) current += nperslice;
2048  //--- copy extra candidates
2049  extra[indextra] = nleft;
2050  extra[indextra+1] = nright;
2051  if (nleft) memcpy(&extra[indextra+2], extra_left, nleft*sizeof(Int_t));
2052  if (nright) memcpy(&extra[indextra+2+nleft], extra_right, nright*sizeof(Int_t));
2053  indextra += 2+nleft+nright;
2054  }
2055  if (fIndcZ) delete [] fIndcZ;
2056  fNz = current;
2057  fIndcZ = new UChar_t[current];
2058  memcpy(fIndcZ, &ind[0], current*sizeof(UChar_t));
2059  if (fExtraZ) delete [] fExtraZ;
2060  fNez = indextra;
2061  if (indextra>nmaxslices*4) printf("Woops!!!\n");
2062  fExtraZ = new Int_t[indextra];
2063  memcpy(fExtraZ, extra, indextra*sizeof(Int_t));
2064  }
2065  delete [] boundaries;
2066  delete [] index;
2067  delete [] temp;
2068  delete [] ind;
2069  delete [] extra;
2070  delete [] extra_left;
2071  delete [] extra_right;
2072 
2073  if ((!fPriority[0]) && (!fPriority[1]) && (!fPriority[2])) {
2074  SetInvalid();
2075  if (nd>1) Error("SortAll", "Volume %s: Cannot make slices on any axis", fVolume->GetName());
2076  }
2077 }
2078 
2079 ////////////////////////////////////////////////////////////////////////////////
2080 /// Print the voxels.
2081 
2083 {
2084  if (NeedRebuild()) {
2085  TGeoVoxelFinder *vox = (TGeoVoxelFinder*)this;
2086  vox->Voxelize();
2087  fVolume->FindOverlaps();
2088  }
2089  Int_t id, i;
2090  Int_t nd = fVolume->GetNdaughters();
2091  printf("Voxels for volume %s (nd=%i)\n", fVolume->GetName(), fVolume->GetNdaughters());
2092  printf("priority : x=%i y=%i z=%i\n", fPriority[0], fPriority[1], fPriority[2]);
2093 // return;
2094  Int_t nextra;
2095  Int_t nbytes = 1+((fVolume->GetNdaughters()-1)>>3);
2096  UChar_t byte, bit;
2097  UChar_t *slice;
2098  printf("XXX\n");
2099  if (fPriority[0]==2) {
2100  for (id=0; id<fIbx; id++) {
2101  printf("%15.10f\n",fXb[id]);
2102  if (id == (fIbx-1)) break;
2103  printf("slice %i : %i\n", id, fNsliceX[id]);
2104  if (fNsliceX[id]) {
2105  slice = &fIndcX[fOBx[id]];
2106  for (i=0; i<nbytes; i++) {
2107  byte = slice[i];
2108  for (bit=0; bit<8; bit++) {
2109  if (byte & (1<<bit)) printf(" %i ", 8*i+bit);
2110  }
2111  }
2112  printf("\n");
2113  }
2114  GetExtraX(id,kTRUE,nextra);
2115  printf(" extra_about_left = %i\n", nextra);
2116  GetExtraX(id,kFALSE,nextra);
2117  printf(" extra_about_right = %i\n", nextra);
2118  }
2119  } else if (fPriority[0]==1) {
2120  printf("%15.10f\n",fXb[0]);
2121  for (id=0; id<nd; id++) printf(" %i ",id);
2122  printf("\n");
2123  printf("%15.10f\n",fXb[1]);
2124  }
2125  printf("YYY\n");
2126  if (fPriority[1]==2) {
2127  for (id=0; id<fIby; id++) {
2128  printf("%15.10f\n", fYb[id]);
2129  if (id == (fIby-1)) break;
2130  printf("slice %i : %i\n", id, fNsliceY[id]);
2131  if (fNsliceY[id]) {
2132  slice = &fIndcY[fOBy[id]];
2133  for (i=0; i<nbytes; i++) {
2134  byte = slice[i];
2135  for (bit=0; bit<8; bit++) {
2136  if (byte & (1<<bit)) printf(" %i ", 8*i+bit);
2137  }
2138  }
2139  }
2140  GetExtraY(id,kTRUE,nextra);
2141  printf(" extra_about_left = %i\n", nextra);
2142  GetExtraY(id,kFALSE,nextra);
2143  printf(" extra_about_right = %i\n", nextra);
2144  }
2145  } else if (fPriority[1]==1) {
2146  printf("%15.10f\n",fYb[0]);
2147  for (id=0; id<nd; id++) printf(" %i ",id);
2148  printf("\n");
2149  printf("%15.10f\n",fYb[1]);
2150  }
2151 
2152  printf("ZZZ\n");
2153  if (fPriority[2]==2) {
2154  for (id=0; id<fIbz; id++) {
2155  printf("%15.10f\n", fZb[id]);
2156  if (id == (fIbz-1)) break;
2157  printf("slice %i : %i\n", id, fNsliceZ[id]);
2158  if (fNsliceZ[id]) {
2159  slice = &fIndcZ[fOBz[id]];
2160  for (i=0; i<nbytes; i++) {
2161  byte = slice[i];
2162  for (bit=0; bit<8; bit++) {
2163  if (byte & (1<<bit)) printf(" %i ", 8*i+bit);
2164  }
2165  }
2166  printf("\n");
2167  }
2168  GetExtraZ(id,kTRUE,nextra);
2169  printf(" extra_about_left = %i\n", nextra);
2170  GetExtraZ(id,kFALSE,nextra);
2171  printf(" extra_about_right = %i\n", nextra);
2172  }
2173  } else if (fPriority[2]==1) {
2174  printf("%15.10f\n",fZb[0]);
2175  for (id=0; id<nd; id++) printf(" %i ",id);
2176  printf("\n");
2177  printf("%15.10f\n",fZb[1]);
2178  }
2179 }
2180 
2181 ////////////////////////////////////////////////////////////////////////////////
2182 /// print the voxel containing point
2183 
2185 {
2186  if (NeedRebuild()) {
2187  TGeoVoxelFinder *vox = (TGeoVoxelFinder*)this;
2188  vox->Voxelize();
2189  fVolume->FindOverlaps();
2190  }
2191  Int_t im=0;
2192  if (fPriority[0]) {
2193  im = TMath::BinarySearch(fIbx, fXb, point[0]);
2194  if ((im==-1) || (im==fIbx-1)) {
2195  printf("Voxel X limits: OUT\n");
2196  } else {
2197  printf("Voxel X limits: %g %g\n", fXb[im], fXb[im+1]);
2198  }
2199  }
2200  if (fPriority[1]) {
2201  im = TMath::BinarySearch(fIby, fYb, point[1]);
2202  if ((im==-1) || (im==fIby-1)) {
2203  printf("Voxel Y limits: OUT\n");
2204  } else {
2205  printf("Voxel Y limits: %g %g\n", fYb[im], fYb[im+1]);
2206  }
2207  }
2208  if (fPriority[2]) {
2209  im = TMath::BinarySearch(fIbz, fZb, point[2]);
2210  if ((im==-1) || (im==fIbz-1)) {
2211  printf("Voxel Z limits: OUT\n");
2212  } else {
2213  printf("Voxel Z limits: %g %g\n", fZb[im], fZb[im+1]);
2214  }
2215  }
2216 }
2217 ////////////////////////////////////////////////////////////////////////////////
2218 /// Voxelize attached volume according to option
2219 /// If the volume is an assembly, make sure the bbox is computed.
2220 
2222 {
2224  Int_t nd = fVolume->GetNdaughters();
2225  TGeoVolume *vd;
2226  for (Int_t i=0; i<nd; i++) {
2227  vd = fVolume->GetNode(i)->GetVolume();
2228  if (vd->IsAssembly()) vd->GetShape()->ComputeBBox();
2229  }
2230  BuildVoxelLimits();
2231  SortAll();
2233 }
2234 ////////////////////////////////////////////////////////////////////////////////
2235 /// Stream an object of class TGeoVoxelFinder.
2236 
2237 void TGeoVoxelFinder::Streamer(TBuffer &R__b)
2238 {
2239  if (R__b.IsReading()) {
2240  UInt_t R__s, R__c;
2241  Version_t R__v = R__b.ReadVersion(&R__s, &R__c);
2242  if (R__v > 2) {
2243  R__b.ReadClassBuffer(TGeoVoxelFinder::Class(), this, R__v, R__s, R__c);
2244  return;
2245  }
2246  // Process old versions of the voxel finder. Just read the data
2247  // from the buffer in a temp variable then mark voxels as garbage.
2248  UChar_t *dummy = new UChar_t[R__c-2];
2249  R__b.ReadFastArray(dummy, R__c-2);
2250  delete [] dummy;
2251  SetInvalid(kTRUE);
2252  } else {
2254  }
2255 }
TGeoVoxelFinder::FindOverlaps
virtual void FindOverlaps(Int_t inode) const
create the list of nodes for which the bboxes overlap with inode's bbox
Definition: TGeoVoxelFinder.cxx:296
TGeoVoxelFinder::PrintVoxelLimits
void PrintVoxelLimits(const Double_t *point) const
print the voxel containing point
Definition: TGeoVoxelFinder.cxx:2184
TGeoVoxelFinder::fNoy
Int_t fNoy
Definition: TGeoVoxelFinder.h:40
TGeoVoxelFinder::fNz
Int_t fNz
Definition: TGeoVoxelFinder.h:47
ymax
float ymax
Definition: THbookFile.cxx:95
kTRUE
const Bool_t kTRUE
Definition: RtypesCore.h:100
TGeoVoxelFinder::GetCheckList
virtual Int_t * GetCheckList(const Double_t *point, Int_t &nelem, TGeoStateInfo &td)
get the list of daughter indices for which point is inside their bbox
Definition: TGeoVoxelFinder.cxx:1205
Version_t
short Version_t
Definition: RtypesCore.h:65
TGeoVolume::GetShape
TGeoShape * GetShape() const
Definition: TGeoVolume.h:188
TGeoVolume::GetNdaughters
Int_t GetNdaughters() const
Definition: TGeoVolume.h:349
Option_t
const char Option_t
Definition: RtypesCore.h:66
TGeoVoxelFinder::fExtraX
Int_t * fExtraX
Definition: TGeoVoxelFinder.h:59
gGeoManager
R__EXTERN TGeoManager * gGeoManager
Definition: TGeoManager.h:602
TGeoVoxelFinder::IsSafeVoxel
Bool_t IsSafeVoxel(const Double_t *point, Int_t inode, Double_t minsafe) const
Computes squared distance from POINT to the voxel(s) containing node INODE.
Definition: TGeoVoxelFinder.cxx:227
TGeoVoxelFinder::fIbz
Int_t fIbz
Definition: TGeoVoxelFinder.h:37
TBuffer::ReadFastArray
virtual void ReadFastArray(Bool_t *b, Int_t n)=0
ClassImp
#define ClassImp(name)
Definition: Rtypes.h:364
TGeoVoxelFinder::Efficiency
virtual Double_t Efficiency()
Compute voxelization efficiency.
Definition: TGeoVoxelFinder.cxx:249
TGeoVoxelFinder::DaughterToMother
void DaughterToMother(Int_t id, const Double_t *local, Double_t *master) const
convert a point from the local reference system of node id to reference system of mother volume
Definition: TGeoVoxelFinder.cxx:217
TBuffer::ReadClassBuffer
virtual Int_t ReadClassBuffer(const TClass *cl, void *pointer, const TClass *onfile_class=0)=0
xmax
float xmax
Definition: THbookFile.cxx:95
TGeoVoxelFinder::SortAll
void SortAll(Option_t *option="")
order bounding boxes along x, y, z
Definition: TGeoVoxelFinder.cxx:1658
TObject::Error
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:893
TGeoNode::GetVolume
TGeoVolume * GetVolume() const
Definition: TGeoNode.h:97
TGeoStateInfo::fVoxCurrent
Int_t fVoxCurrent
Definition: TGeoStateInfo.h:34
TGeoVoxelFinder::fNboxes
Int_t fNboxes
Definition: TGeoVoxelFinder.h:38
TObject::Fatal
virtual void Fatal(const char *method, const char *msgfmt,...) const
Issue fatal error message.
Definition: TObject.cxx:921
TGeoVoxelFinder::GetNextVoxel
virtual Int_t * GetNextVoxel(const Double_t *point, const Double_t *dir, Int_t &ncheck, TGeoStateInfo &td)
get the list of new candidates for the next voxel crossed by current ray printf("### GetNextVoxel\n")...
Definition: TGeoVoxelFinder.cxx:1366
box
void box(Int_t pat, Double_t x1, Double_t y1, Double_t x2, Double_t y2)
Definition: fillpatterns.C:1
TGeoVoxelFinder::GetVoxelCandidates
Int_t * GetVoxelCandidates(Int_t i, Int_t j, Int_t k, Int_t &ncheck, TGeoStateInfo &td)
get the list of candidates in voxel (i,j,k) - no check
Definition: TGeoVoxelFinder.cxx:1302
TGeoVoxelFinder::GetIndices
Bool_t GetIndices(const Double_t *point, TGeoStateInfo &td)
Get indices for current slices on x, y, z.
Definition: TGeoVoxelFinder.cxx:348
TGeoVoxelFinder::TGeoVoxelFinder
TGeoVoxelFinder()
Default constructor.
Definition: TGeoVoxelFinder.cxx:39
TGeoVoxelFinder::fIby
Int_t fIby
Definition: TGeoVoxelFinder.h:36
TMath::Abs
Short_t Abs(Short_t d)
Definition: TMathBase.h:120
TBuffer
Buffer base class used for serializing objects.
Definition: TBuffer.h:43
TGeoVoxelFinder::BuildVoxelLimits
void BuildVoxelLimits()
build the array of bounding boxes of the nodes inside
Definition: TGeoVoxelFinder.cxx:172
TMath::Sort
void Sort(Index n, const Element *a, Index *index, Bool_t down=kTRUE)
Definition: TMathBase.h:362
TGeoVoxelFinder::Union
Bool_t Union(Int_t n1, UChar_t *array1, TGeoStateInfo &td)
make union of older bits with new array printf("Union - one slice\n");
Definition: TGeoVoxelFinder.cxx:1460
TGeoVoxelFinder::fOBx
Int_t * fOBx
Definition: TGeoVoxelFinder.h:53
TGeoBBox.h
TGeoNode::GetMatrix
virtual TGeoMatrix * GetMatrix() const =0
TGeoNode
A node represent a volume positioned inside another.They store links to both volumes and to the TGeoM...
Definition: TGeoNode.h:41
TGeoVoxelFinder::fOEx
Int_t * fOEx
Definition: TGeoVoxelFinder.h:56
bool
TGeoStateInfo::fVoxInvdir
Double_t fVoxInvdir[3]
Definition: TGeoStateInfo.h:39
id
XFontStruct * id
Definition: TGX11.cxx:109
TGeoVoxelFinder::fVolume
TGeoVolume * fVolume
Definition: TGeoVoxelFinder.h:33
TGeoVoxelFinder::GetNextCandidates
virtual Int_t * GetNextCandidates(const Double_t *point, Int_t &ncheck, TGeoStateInfo &td)
Returns list of new candidates in next voxel.
Definition: TGeoVoxelFinder.cxx:510
TGeoStateInfo::fVoxCheckList
Int_t * fVoxCheckList
Definition: TGeoStateInfo.h:35
TGeoVoxelFinder::fBoxes
Double_t * fBoxes
Definition: TGeoVoxelFinder.h:49
TGeoVoxelFinder::SetInvalid
void SetInvalid(Bool_t flag=kTRUE)
Definition: TGeoVoxelFinder.h:110
TGeoVoxelFinder::fNex
Int_t fNex
Definition: TGeoVoxelFinder.h:42
TGeoStateInfo::fVoxSlices
Int_t fVoxSlices[3]
Definition: TGeoStateInfo.h:37
TGeoStateInfo
Statefull info for the current geometry level.
Definition: TGeoStateInfo.h:21
TMath::LocMin
Long64_t LocMin(Long64_t n, const T *a)
Return index of array with the minimum element.
Definition: TMath.h:972
TGeoVoxelFinder::fXb
Double_t * fXb
Definition: TGeoVoxelFinder.h:50
TGeoVoxelFinder::fNsliceX
Int_t * fNsliceX
Definition: TGeoVoxelFinder.h:62
TBuffer.h
TGeoVoxelFinder::fExtraY
Int_t * fExtraY
Definition: TGeoVoxelFinder.h:60
TGeoVoxelFinder::Print
virtual void Print(Option_t *option="") const
Print the voxels.
Definition: TGeoVoxelFinder.cxx:2082
xmin
float xmin
Definition: THbookFile.cxx:95
TGeoShape::ComputeBBox
virtual void ComputeBBox()=0
TGeoVoxelFinder::SortCrossedVoxels
virtual void SortCrossedVoxels(const Double_t *point, const Double_t *dir, TGeoStateInfo &td)
get the list in the next voxel crossed by a ray
Definition: TGeoVoxelFinder.cxx:1071
TBuffer::WriteClassBuffer
virtual Int_t WriteClassBuffer(const TClass *cl, void *pointer)=0
TGeoVoxelFinder::fOEz
Int_t * fOEz
Definition: TGeoVoxelFinder.h:58
kFALSE
const Bool_t kFALSE
Definition: RtypesCore.h:101
TGeoVoxelFinder::fOEy
Int_t * fOEy
Definition: TGeoVoxelFinder.h:57
TGeoVoxelFinder::fYb
Double_t * fYb
Definition: TGeoVoxelFinder.h:51
TGeoVoxelFinder::fNx
Int_t fNx
Definition: TGeoVoxelFinder.h:45
TGeoVolume::SetCylVoxels
void SetCylVoxels(Bool_t flag=kTRUE)
Definition: TGeoVolume.h:213
TGeoVoxelFinder::GetExtraZ
Int_t * GetExtraZ(Int_t islice, Bool_t left, Int_t &nextra) const
Return the list of extra candidates in a given Z slice compared to another (left or right)
Definition: TGeoVoxelFinder.cxx:431
UInt_t
unsigned int UInt_t
Definition: RtypesCore.h:46
TGeoVoxelFinder::fOBz
Int_t * fOBz
Definition: TGeoVoxelFinder.h:55
TGeoBBox
Box class.
Definition: TGeoBBox.h:18
TGeoStateInfo::fVoxNcandidates
Int_t fVoxNcandidates
Definition: TGeoStateInfo.h:33
TGeoVoxelFinder::SetNeedRebuild
void SetNeedRebuild(Bool_t flag=kTRUE)
Definition: TGeoVoxelFinder.h:111
TGeoVoxelFinder::Voxelize
virtual void Voxelize(Option_t *option="")
Voxelize attached volume according to option If the volume is an assembly, make sure the bbox is comp...
Definition: TGeoVoxelFinder.cxx:2221
TGeoMatrix::LocalToMaster
virtual void LocalToMaster(const Double_t *local, Double_t *master) const
convert a point by multiplying its column vector (x, y, z, 1) to matrix inverse
Definition: TGeoMatrix.cxx:339
TGeoVolume::IsAssembly
virtual Bool_t IsAssembly() const
Returns true if the volume is an assembly or a scaled assembly.
Definition: TGeoVolume.cxx:1705
TGeoVoxelFinder.h
TBuffer::ReadVersion
virtual Version_t ReadVersion(UInt_t *start=0, UInt_t *bcnt=0, const TClass *cl=0)=0
TGeoVolume::FindOverlaps
void FindOverlaps() const
loop all nodes marked as overlaps and find overlapping brothers
Definition: TGeoVolume.cxx:2090
TGeoVoxelFinder::fNsliceY
Int_t * fNsliceY
Definition: TGeoVoxelFinder.h:63
byte
unsigned char byte
Definition: gifdecode.c:10
TMath::Min
Short_t Min(Short_t a, Short_t b)
Definition: TMathBase.h:180
unsigned int
TGeoMatrix.h
TGeoVoxelFinder::fNoz
Int_t fNoz
Definition: TGeoVoxelFinder.h:41
ymin
float ymin
Definition: THbookFile.cxx:95
TMath::BinarySearch
Long64_t BinarySearch(Long64_t n, const T *array, T value)
Definition: TMathBase.h:278
TGeoVoxelFinder::fIndcX
UChar_t * fIndcX
Definition: TGeoVoxelFinder.h:65
TGeoManager.h
TBuffer::IsReading
Bool_t IsReading() const
Definition: TBuffer.h:86
TGeoVoxelFinder::fOBy
Int_t * fOBy
Definition: TGeoVoxelFinder.h:54
TGeoStateInfo::fVoxBits1
UChar_t * fVoxBits1
Definition: TGeoStateInfo.h:36
TGeoVoxelFinder::fIndcZ
UChar_t * fIndcZ
Definition: TGeoVoxelFinder.h:67
TGeoVoxelFinder::fNox
Int_t fNox
Definition: TGeoVoxelFinder.h:39
Double_t
double Double_t
Definition: RtypesCore.h:59
TGeoMatrix
Geometrical transformation package.
Definition: TGeoMatrix.h:41
TGeoVoxelFinder::NeedRebuild
Bool_t NeedRebuild() const
Definition: TGeoVoxelFinder.h:105
TGeoVoxelFinder::~TGeoVoxelFinder
virtual ~TGeoVoxelFinder()
Destructor.
Definition: TGeoVoxelFinder.cxx:127
TGeoVolume::GetNode
TGeoNode * GetNode(const char *name) const
get the pointer to a daughter node
Definition: TGeoVolume.cxx:2060
TGeoVoxelFinder::GetValidExtra
Int_t * GetValidExtra(Int_t *list, Int_t &ncheck, TGeoStateInfo &td)
Get extra candidates that are not contained in current check list.
Definition: TGeoVoxelFinder.cxx:449
UChar_t
unsigned char UChar_t
Definition: RtypesCore.h:38
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:326
TGeoVoxelFinder::fExtraZ
Int_t * fExtraZ
Definition: TGeoVoxelFinder.h:61
TGeoVoxelFinder::Intersect
Bool_t Intersect(Int_t n1, UChar_t *array1, Int_t &nf, Int_t *result)
return the list of nodes corresponding to one array of bits
Definition: TGeoVoxelFinder.cxx:1389
TGeoVoxelFinder::fZb
Double_t * fZb
Definition: TGeoVoxelFinder.h:52
TGeoNode::SetOverlaps
void SetOverlaps(Int_t *ovlp, Int_t novlp)
set the list of overlaps for this node (ovlp must be created with operator new)
Definition: TGeoNode.cxx:653
TGeoStateInfo::fVoxLimits
Double_t fVoxLimits[3]
Definition: TGeoStateInfo.h:40
TGeoShape::Tolerance
static Double_t Tolerance()
Definition: TGeoShape.h:91
TGeoVoxelFinder::fNey
Int_t fNey
Definition: TGeoVoxelFinder.h:43
TGeoVoxelFinder::fNsliceZ
Int_t * fNsliceZ
Definition: TGeoVoxelFinder.h:64
TGeoVoxelFinder::IntersectAndStore
Bool_t IntersectAndStore(Int_t n1, UChar_t *array1, TGeoStateInfo &td)
return the list of nodes corresponding to one array of bits
Definition: TGeoVoxelFinder.cxx:1418
TGeoVoxelFinder::GetExtraY
Int_t * GetExtraY(Int_t islice, Bool_t left, Int_t &nextra) const
Return the list of extra candidates in a given Y slice compared to another (left or right)
Definition: TGeoVoxelFinder.cxx:412
TNamed::GetName
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
TGeoStateInfo.h
TGeoShape::Big
static Double_t Big()
Definition: TGeoShape.h:88
TGeoStateInfo::fVoxInc
Int_t fVoxInc[3]
Definition: TGeoStateInfo.h:38
pt
TPaveText * pt
Definition: entrylist_figure1.C:7
Class
void Class()
Definition: Class.C:29
TGeoVoxelFinder::fIbx
Int_t fIbx
Definition: TGeoVoxelFinder.h:35
TGeoVoxelFinder::fPriority
Int_t fPriority[3]
Definition: TGeoVoxelFinder.h:48
TGeoVoxelFinder::GetNcandidates
Int_t GetNcandidates(TGeoStateInfo &td) const
Definition: TGeoVoxelFinder.cxx:156
TGeoVoxelFinder
Finder class handling voxels.
Definition: TGeoVoxelFinder.h:21
TMath::E
constexpr Double_t E()
Base of natural log:
Definition: TMath.h:96
TGeoVolume
TGeoVolume, TGeoVolumeMulti, TGeoVolumeAssembly are the volume classes.
Definition: TGeoVolume.h:49
TGeoNode.h
TGeoVoxelFinder::GetExtraX
Int_t * GetExtraX(Int_t islice, Bool_t left, Int_t &nextra) const
Return the list of extra candidates in a given X slice compared to another (left or right)
Definition: TGeoVoxelFinder.cxx:393
TGeoVoxelFinder::fIndcY
UChar_t * fIndcY
Definition: TGeoVoxelFinder.h:66
TMath.h
TGeoVoxelFinder::fNy
Int_t fNy
Definition: TGeoVoxelFinder.h:46
int
TGeoManager::GetStep
Double_t GetStep() const
Definition: TGeoManager.h:408
TGeoVoxelFinder::fNez
Int_t fNez
Definition: TGeoVoxelFinder.h:44