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