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