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