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