Logo ROOT  

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