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