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