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