Logo ROOT  

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