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