    // @(#)root/geom:$Name:  $:$Id: TGeoVoxelFinder.cxx,v 1.24 2005/06/13 12:17:32 brun Exp $
    // Author: Andrei Gheata   04/02/02

    /*************************************************************************
     * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers.               *
     * All rights reserved.                                                  *
     *                                                                       *
     * For the licensing terms see $ROOTSYS/LICENSE.                         *
     * For the list of contributors see $ROOTSYS/README/CREDITS.             *
     *************************************************************************/

    ////////////////////////////////////////////////////////////////////////////////
    // Full description with examples and pictures
    //
    //
    //
    //
    //

/* */ //

#include "TObject.h" #include "TGeoMatrix.h" #include "TGeoBBox.h" #include "TGeoNode.h" #include "TGeoManager.h" #include "TGeoVoxelFinder.h" /************************************************************************* * TGeoVoxelFinder - finder class handling voxels * *************************************************************************/ ClassImp(TGeoVoxelFinder) //----------------------------------------------------------------------------- TGeoVoxelFinder::TGeoVoxelFinder() { // Default constructor fVolume = 0; fNboxes = 0; fIbx = 0; fIby = 0; fIbz = 0; fNox = 0; fNoy = 0; fNoz = 0; fNex = 0; fNey = 0; fNez = 0; fNx = 0; fNy = 0; fNz = 0; fBoxes = 0; fXb = 0; fYb = 0; fZb = 0; fOBx = 0; fOBy = 0; fOBz = 0; fOEx = 0; fOEy = 0; fOEz = 0; fIndX = 0; fIndY = 0; fIndZ = 0; fExtraX = 0; fExtraY = 0; fExtraZ = 0; fCheckList = 0; fNcandidates = 0; fCurrentVoxel = 0; fBits1 = 0; SetInvalid(kFALSE); } //----------------------------------------------------------------------------- TGeoVoxelFinder::TGeoVoxelFinder(TGeoVolume *vol) { // Default constructor if (!vol) return; fVolume = vol; fVolume->SetCylVoxels(kFALSE); fNboxes = 0; fIbx = 0; fIby = 0; fIbz = 0; fNox = 0; fNoy = 0; fNoz = 0; fNex = 0; fNey = 0; fNez = 0; fNx = 0; fNy = 0; fNz = 0; fBoxes = 0; fXb = 0; fYb = 0; fZb = 0; fOBx = 0; fOBy = 0; fOBz = 0; fOEx = 0; fOEy = 0; fOEz = 0; fIndX = 0; fIndY = 0; fIndZ = 0; fExtraX = 0; fExtraY = 0; fExtraZ = 0; fCheckList = 0; fNcandidates = 0; fCurrentVoxel = 0; fBits1 = 0; } //----------------------------------------------------------------------------- TGeoVoxelFinder::~TGeoVoxelFinder() { // Destructor // printf("deleting finder of %s\n", fVolume->GetName()); if (fOBx) delete [] fOBx; if (fOBy) delete [] fOBy; if (fOBz) delete [] fOBz; // printf("OBx OBy OBz...\n"); if (fBoxes) delete [] fBoxes; // printf("boxes...\n"); if (fXb) delete [] fXb; if (fYb) delete [] fYb; if (fZb) delete [] fZb; // printf("Xb Yb Zb...\n"); if (fIndX) delete [] fIndX; if (fIndY) delete [] fIndY; if (fIndZ) delete [] fIndZ; if (fExtraX) delete [] fExtraX; if (fExtraY) delete [] fExtraY; if (fExtraZ) delete [] fExtraZ; // printf("IndX IndY IndZ...\n"); if (fCheckList) delete [] fCheckList; if (fBits1) delete [] fBits1; // printf("checklist...\n"); } //----------------------------------------------------------------------------- void TGeoVoxelFinder::BuildVoxelLimits() { // build the array of bounding boxes of the nodes inside Int_t nd = fVolume->GetNdaughters(); if (!nd) return; //printf("building boxes for %s nd=%i\n", fVolume->GetName(), nd); Int_t id; TGeoNode *node; if (fBoxes) delete [] fBoxes; if (fBits1) delete [] fBits1; fBits1 = new UChar_t[1+((nd-1)>>3)]; fNboxes = 6*nd; fBoxes = new Double_t[fNboxes]; if (fCheckList) delete [] fCheckList; fCheckList = new Int_t[nd]; Double_t vert[24]; Double_t pt[3]; Double_t xyz[6]; // printf("boundaries for %s :\n", GetName()); TGeoBBox *box = 0; for (id=0; id<nd; id++) { node = fVolume->GetNode(id); // if (!strcmp(node->ClassName(), "TGeoNodeOffset") continue; box = (TGeoBBox*)node->GetVolume()->GetShape(); box->SetBoxPoints(&vert[0]); for (Int_t point=0; point<8; point++) { DaughterToMother(id, &vert[3*point], &pt[0]); if (!point) { xyz[0] = xyz[1] = pt[0]; xyz[2] = xyz[3] = pt[1]; xyz[4] = xyz[5] = pt[2]; continue; } for (Int_t j=0; j<3; j++) { if (pt[j] < xyz[2*j]) xyz[2*j]=pt[j]; if (pt[j] > xyz[2*j+1]) xyz[2*j+1]=pt[j]; } } fBoxes[6*id] = 0.5*(xyz[1]-xyz[0]); // dX fBoxes[6*id+1] = 0.5*(xyz[3]-xyz[2]); // dY fBoxes[6*id+2] = 0.5*(xyz[5]-xyz[4]); // dZ fBoxes[6*id+3] = 0.5*(xyz[0]+xyz[1]); // Ox fBoxes[6*id+4] = 0.5*(xyz[2]+xyz[3]); // Oy fBoxes[6*id+5] = 0.5*(xyz[4]+xyz[5]); // Oz } } //----------------------------------------------------------------------------- void TGeoVoxelFinder::CreateCheckList() { // Initializes check list. Int_t nd = fVolume->GetNdaughters(); if (!nd) return; if (!fCheckList) fCheckList = new Int_t[nd]; if (!fBits1) fBits1 = new UChar_t[1+((nd-1)>>3)]; } //----------------------------------------------------------------------------- void TGeoVoxelFinder::DaughterToMother(Int_t id, Double_t *local, Double_t *master) const { // convert a point from the local reference system of node id to reference // system of mother volume TGeoMatrix *mat = fVolume->GetNode(id)->GetMatrix(); mat->LocalToMaster(local, master); } //----------------------------------------------------------------------------- Bool_t TGeoVoxelFinder::IsSafeVoxel(Double_t *point, Int_t inode, Double_t minsafe) const { // Computes squared distance from POINT to the voxel(s) containing node INODE. Returns 0 // if POINT inside voxel(s). Double_t dxyz, minsafe2=minsafe*minsafe; Int_t ist = 6*inode; Int_t i; Double_t rsq = 0; for (i=0; i<3; i++) { dxyz = TMath::Abs(point[i]-fBoxes[ist+i+3])-fBoxes[ist+i]; if (dxyz>-1E-6) rsq+=dxyz*dxyz; if (rsq >= minsafe2) return kTRUE; } return kFALSE; } //----------------------------------------------------------------------------- Double_t TGeoVoxelFinder::Efficiency() { //--- Compute voxelization efficiency. printf("Voxelization efficiency for %s\n", fVolume->GetName()); Double_t nd = Double_t(fVolume->GetNdaughters()); Double_t eff = 0; Double_t effslice = 0; Int_t id; if (fPriority[0]) { for (id=0; id<fIbx-1; id++) { // loop on boundaries effslice += fIndX[fOBx[id]]; } if (effslice != 0) effslice = nd/effslice; else printf("Woops : slice X\n"); } printf("X efficiency : %g\n", effslice); eff += effslice; effslice = 0; if (fPriority[1]) { for (id=0; id<fIby-1; id++) { // loop on boundaries effslice += fIndY[fOBy[id]]; } if (effslice != 0) effslice = nd/effslice; else printf("Woops : slice X\n"); } printf("Y efficiency : %g\n", effslice); eff += effslice; effslice = 0; if (fPriority[2]) { for (id=0; id<fIbz-1; id++) { // loop on boundaries effslice += fIndZ[fOBz[id]]; } if (effslice != 0) effslice = nd/effslice; else printf("Woops : slice X\n"); } printf("Z efficiency : %g\n", effslice); eff += effslice; eff /= 3.; printf("Total efficiency : %g\n", eff); return eff; } //----------------------------------------------------------------------------- void TGeoVoxelFinder::FindOverlaps(Int_t inode) const { // create the list of nodes for which the bboxes overlap with inode's bbox if (!fBoxes) return; Double_t xmin, xmax, ymin, ymax, zmin, zmax; Double_t xmin1, xmax1, ymin1, ymax1, zmin1, zmax1; Double_t ddx1, ddx2; Int_t nd = fVolume->GetNdaughters(); Int_t *ovlps = 0; Int_t *otmp = new Int_t[nd-1]; Int_t novlp = 0; TGeoNode *node = fVolume->GetNode(inode); // printf("Finding overlaps for %s\n", node->GetName()); xmin = fBoxes[6*inode+3] - fBoxes[6*inode]; xmax = fBoxes[6*inode+3] + fBoxes[6*inode]; ymin = fBoxes[6*inode+4] - fBoxes[6*inode+1]; ymax = fBoxes[6*inode+4] + fBoxes[6*inode+1]; zmin = fBoxes[6*inode+5] - fBoxes[6*inode+2]; zmax = fBoxes[6*inode+5] + fBoxes[6*inode+2]; // printf("overlaps for MANY node %s\n", node->GetName()); // printf("xmin=%g xmax=%g\n", xmin, xmax); // printf("ymin=%g ymax=%g\n", ymin, ymax); // printf("zmin=%g zmax=%g\n", zmin, zmax); //TGeoNode *node1; // loop on brothers for (Int_t ib=0; ib<nd; ib++) { if (ib == inode) continue; // everyone overlaps with itself //node1 = fVolume->GetNode(ib); xmin1 = fBoxes[6*ib+3] - fBoxes[6*ib]; xmax1 = fBoxes[6*ib+3] + fBoxes[6*ib]; ymin1 = fBoxes[6*ib+4] - fBoxes[6*ib+1]; ymax1 = fBoxes[6*ib+4] + fBoxes[6*ib+1]; zmin1 = fBoxes[6*ib+5] - fBoxes[6*ib+2]; zmax1 = fBoxes[6*ib+5] + fBoxes[6*ib+2]; // printf(" node %s\n", node1->GetName()); // printf(" xmin1=%g xmax1=%g\n", xmin1, xmax1); // printf(" ymin1=%g ymax1=%g\n", ymin1, ymax1); // printf(" zmin1=%g zmax1=%g\n", zmin1, zmax1); ddx1 = xmin1-xmax; ddx2 = xmax1-xmin; if ((ddx1>-1E-10)||(ddx2<1E-10)) continue; ddx1 = ymin1-ymax; ddx2 = ymax1-ymin; if ((ddx1>-1E-10)||(ddx2<1E-10)) continue; ddx1 = zmin1-zmax; ddx2 = zmax1-zmin; if ((ddx1>-1E-10)||(ddx2<1E-10)) continue; otmp[novlp++] = ib; } if (!novlp) { // printf("---no overlaps for MANY node %s\n", node->GetName()); node->SetOverlaps(ovlps, 0); return; } ovlps = new Int_t[novlp]; memcpy(ovlps, otmp, novlp*sizeof(Int_t)); delete [] otmp; node->SetOverlaps(ovlps, novlp); // printf("Overlaps for MANY node %s : %i\n", node->GetName(), novlp); } //----------------------------------------------------------------------------- Bool_t TGeoVoxelFinder::GetIndices(Double_t *point) { // Getindices for current slices on x, y, z fSlices[0] = -2; // -2 means 'all daughters in slice' fSlices[1] = -2; fSlices[2] = -2; Bool_t flag=kTRUE; if (fPriority[0]) { fSlices[0] = TMath::BinarySearch(fIbx, fXb, point[0]); if ((fSlices[0]<0) || (fSlices[0]>=fIbx-1)) { // outside slices flag=kFALSE; } else { if (fPriority[0]==2) { // nothing in current slice if (!fIndX[fOBx[fSlices[0]]]) flag = kFALSE; } } } if (fPriority[1]) { fSlices[1] = TMath::BinarySearch(fIby, fYb, point[1]); if ((fSlices[1]<0) || (fSlices[1]>=fIby-1)) { // outside slices flag=kFALSE; } else { if (fPriority[1]==2) { // nothing in current slice if (!fIndY[fOBy[fSlices[1]]]) flag = kFALSE; } } } if (fPriority[2]) { fSlices[2] = TMath::BinarySearch(fIbz, fZb, point[2]); if ((fSlices[2]<0) || (fSlices[2]>=fIbz-1)) return kFALSE; if (fPriority[2]==2) { // nothing in current slice if (!fIndZ[fOBz[fSlices[2]]]) return kFALSE; } } return flag; } //----------------------------------------------------------------------------- Int_t *TGeoVoxelFinder::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) Int_t *list = 0; nextra = 0; if (fPriority[0]!=2) return list; if (left) { nextra = fExtraX[fOEx[islice]]; list = &fExtraX[fOEx[islice]+2]; } else { nextra = fExtraX[fOEx[islice]+1]; list = &fExtraX[fOEx[islice]+2+fExtraX[fOEx[islice]]]; } return list; } //----------------------------------------------------------------------------- Int_t *TGeoVoxelFinder::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) Int_t *list = 0; nextra = 0; if (fPriority[1]!=2) return list; if (left) { nextra = fExtraY[fOEy[islice]]; list = &fExtraY[fOEy[islice]+2]; } else { nextra = fExtraY[fOEy[islice]+1]; list = &fExtraY[fOEy[islice]+2+fExtraY[fOEy[islice]]]; } return list; } //----------------------------------------------------------------------------- Int_t *TGeoVoxelFinder::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) Int_t *list = 0; nextra = 0; if (fPriority[2]!=2) return list; if (left) { nextra = fExtraZ[fOEz[islice]]; list = &fExtraZ[fOEz[islice]+2]; } else { nextra = fExtraZ[fOEz[islice]+1]; list = &fExtraZ[fOEz[islice]+2+fExtraZ[fOEz[islice]]]; } return list; } //----------------------------------------------------------------------------- Int_t *TGeoVoxelFinder::GetValidExtra(Int_t *list, Int_t &ncheck) { // Get extra candidates that are not contained in current check list // UChar_t *bits = gGeoManager->GetBits(); fNcandidates = 0; Int_t icand; UInt_t bitnumber, loc; UChar_t bit, byte; for (icand=0; icand<ncheck; icand++) { bitnumber = (UInt_t)list[icand]; loc = bitnumber>>3; bit = bitnumber%8; byte = (~fBits1[loc]) & (1<<bit); if (byte) fCheckList[fNcandidates++]=list[icand]; } ncheck = fNcandidates; return fCheckList; } //----------------------------------------------------------------------------- Int_t *TGeoVoxelFinder::GetValidExtra(Int_t /*n1*/, UChar_t *array1, Int_t *list, Int_t &ncheck) { // Get extra candidates that are contained in array1 but not in current check list // UChar_t *bits = gGeoManager->GetBits(); fNcandidates = 0; Int_t icand; UInt_t bitnumber, loc; UChar_t bit, byte; for (icand=0; icand<ncheck; icand++) { bitnumber = (UInt_t)list[icand]; loc = bitnumber>>3; bit = bitnumber%8; byte = (~fBits1[loc]) & array1[loc] & (1<<bit); if (byte) fCheckList[fNcandidates++]=list[icand]; } ncheck = fNcandidates; return fCheckList; } //----------------------------------------------------------------------------- Int_t *TGeoVoxelFinder::GetValidExtra(Int_t /*n1*/, UChar_t *array1, Int_t /*n2*/, UChar_t *array2, Int_t *list, Int_t &ncheck) { // Get extra candidates that are contained in array1 but not in current check list // UChar_t *bits = gGeoManager->GetBits(); fNcandidates = 0; Int_t icand; UInt_t bitnumber, loc; UChar_t bit, byte; for (icand=0; icand<ncheck; icand++) { bitnumber = (UInt_t)list[icand]; loc = bitnumber>>3; bit = bitnumber%8; byte = (~fBits1[loc]) & array1[loc] & array2[loc] & (1<<bit); if (byte) fCheckList[fNcandidates++]=list[icand]; } ncheck = fNcandidates; return fCheckList; } //----------------------------------------------------------------------------- Int_t *TGeoVoxelFinder::GetNextCandidates(Double_t *point, Int_t &ncheck) { // Returns list of new candidates in next voxel. If NULL, nowhere to // go next. ncheck = 0; if (fLimits[0]<0) return 0; if (fLimits[1]<0) return 0; if (fLimits[2]<0) return 0; Int_t dind[3]; // new slices //---> start from old slices memcpy(&dind[0], &fSlices[0], 3*sizeof(Int_t)); Double_t dmin[3]; // distances to get to next X,Y, Z slices. dmin[0] = dmin[1] = dmin[2] = TGeoShape::Big(); //---> max. possible step to be considered Double_t maxstep = TMath::Min(gGeoManager->GetStep(), fLimits[TMath::LocMin(3, fLimits)]); // printf("1- maxstep=%g\n", maxstep); Bool_t isXlimit=kFALSE, isYlimit=kFALSE, isZlimit=kFALSE; Bool_t isForcedX=kFALSE, isForcedY=kFALSE, isForcedZ=kFALSE; Double_t dforced[3]; dforced[0] = dforced[1] = dforced[2] = TGeoShape::Big(); Int_t iforced = 0; // //---> work on X if (fPriority[0] && fInc[0]) { //---> increment/decrement slice dind[0] += fInc[0]; if (dind[0]<-1) return 0; // outside range if (dind[0]>fIbx-1) return 0; // outside range if (fInc[0]==1) { dmin[0] = (fXb[dind[0]]-point[0])*fInvdir[0]; } else { dmin[0] = (fXb[fSlices[0]]-point[0])*fInvdir[0]; } isXlimit = (dmin[0]>maxstep)?kTRUE:kFALSE; // printf("---> X : priority=%i, slice=%i/%i inc=%i\n", // fPriority[0], fSlices[0], fIbx-2, fInc[0]); // printf("2- step to next X (%i) = %g\n", (Int_t)isXlimit, dmin[0]); //---> check if propagation to next slice on this axis is forced if ((fSlices[0]==-1) || (fSlices[0]==fIbx-1)) { isForcedX = kTRUE; dforced[0] = dmin[0]; iforced++; // printf(" FORCED 1\n"); if (isXlimit) return 0; } else { if (fPriority[0]==2) { // if no candidates in current slice, force next slice if (fIndX[fOBx[fSlices[0]]]==0) { isForcedX = kTRUE; dforced[0] = dmin[0]; iforced++; // printf(" FORCED 2\n"); if (isXlimit) return 0; } } } } else { // no slices on this axis -> bounding box limit // printf(" No slice on X\n"); dmin[0] = fLimits[0]; isXlimit = kTRUE; } //---> work on Y if (fPriority[1] && fInc[1]) { //---> increment/decrement slice dind[1] += fInc[1]; if (dind[1]<-1) return 0; // outside range if (dind[1]>fIby-1) return 0; // outside range if (fInc[1]==1) { dmin[1] = (fYb[dind[1]]-point[1])*fInvdir[1]; } else { dmin[1] = (fYb[fSlices[1]]-point[1])*fInvdir[1]; } isYlimit = (dmin[1]>maxstep)?kTRUE:kFALSE; // printf("---> Y : priority=%i, slice=%i/%i inc=%i\n", // fPriority[1], fSlices[1], fIby-2, fInc[1]); // printf("3- step to next Y (%i) = %g\n", (Int_t)isYlimit, dmin[1]); //---> check if propagation to next slice on this axis is forced if ((fSlices[1]==-1) || (fSlices[1]==fIby-1)) { isForcedY = kTRUE; dforced[1] = dmin[1]; iforced++; // printf(" FORCED 1\n"); if (isYlimit) return 0; } else { if (fPriority[1]==2) { // if no candidates in current slice, force next slice if (fIndY[fOBy[fSlices[1]]]==0) { isForcedY = kTRUE; dforced[1] = dmin[1]; iforced++; // printf(" FORCED 2\n"); if (isYlimit) return 0; } } } } else { // no slices on this axis -> bounding box limit // printf(" No slice on Y\n"); dmin[1] = fLimits[1]; isYlimit = kTRUE; } //---> work on Z if (fPriority[2] && fInc[2]) { //---> increment/decrement slice dind[2] += fInc[2]; if (dind[2]<-1) return 0; // outside range if (dind[2]>fIbz-1) return 0; // outside range if (fInc[2]==1) { dmin[2] = (fZb[dind[2]]-point[2])*fInvdir[2]; } else { dmin[2] = (fZb[fSlices[2]]-point[2])*fInvdir[2]; } isZlimit = (dmin[2]>maxstep)?kTRUE:kFALSE; // printf("---> Z : priority=%i, slice=%i/%i inc=%i\n", // fPriority[2], fSlices[2], fIbz-2, fInc[2]); // printf("4- step to next Z (%i) = %g\n", (Int_t)isZlimit, dmin[2]); //---> check if propagation to next slice on this axis is forced if ((fSlices[2]==-1) || (fSlices[2]==fIbz-1)) { isForcedZ = kTRUE; dforced[2] = dmin[2]; iforced++; // printf(" FORCED 1\n"); if (isZlimit) return 0; } else { if (fPriority[2]==2) { // if no candidates in current slice, force next slice if (fIndZ[fOBz[fSlices[2]]]==0) { isForcedZ = kTRUE; dforced[2] = dmin[2]; iforced++; // printf(" FORCED 2\n"); if (isZlimit) return 0; } } } } else { // no slices on this axis -> bounding box limit // printf(" No slice on Z\n"); dmin[2] = fLimits[2]; isZlimit = kTRUE; } //---> We are done with checking. See which is the closest slice. // First check if some slice is forced Double_t dslice = 0; Int_t islice = 0; if (iforced) { // some slice is forced if (isForcedX) { // X forced dslice = dforced[0]; islice = 0; if (isForcedY) { // X+Y forced if (dforced[1]>dslice) { dslice = dforced[1]; islice = 1; } if (isForcedZ) { // X+Y+Z forced if (dforced[2]>dslice) { dslice = dforced[2]; islice = 2; } } } else { // X forced if (isForcedZ) { // X+Z forced if (dforced[2]>dslice) { dslice = dforced[2]; islice = 2; } } } } else { if (isForcedY) { // Y forced dslice = dforced[1]; islice = 1; if (isForcedZ) { // Y+Z forced if (dforced[2]>dslice) { dslice = dforced[2]; islice = 2; } } } else { // Z forced dslice = dforced[2]; islice = 2; } } } else { // Nothing forced -> get minimum distance islice = TMath::LocMin(3, dmin); dslice = dmin[islice]; if (dslice>=maxstep) { // printf("DSLICE > MAXSTEP -> EXIT\n"); return 0; } } // printf("5- islicenext=%i DSLICE=%g\n", islice, dslice); Double_t xptnew; Int_t *new_list; // list of new candidates UChar_t *slice1 = 0; UChar_t *slice2 = 0; Int_t ndd[2] = {0,0}; Int_t islices = 0; Bool_t left; switch (islice) { case 0: if (isXlimit) return 0; // increment/decrement X slice fSlices[0]=dind[0]; if (iforced) { // we have to recompute Y and Z slices if (dslice>fLimits[1]) return 0; if (dslice>fLimits[2]) return 0; if ((dslice>dmin[1]) && fInc[1]) { xptnew = point[1]+dslice/fInvdir[1]; // printf(" recomputing Y slice, pos=%g\n", xptnew); while (1) { fSlices[1] += fInc[1]; if (fInc[1]==1) { if (fYb[fSlices[1]+1]>=xptnew) break; } else { if (fYb[fSlices[1]]<= xptnew) break; } } // printf(" %i/%i\n", fSlices[1], fIby-2); } if ((dslice>dmin[2]) && fInc[2]) { xptnew = point[2]+dslice/fInvdir[2]; // printf(" recomputing Z slice, pos=%g\n", xptnew); while (1) { fSlices[2] += fInc[2]; if (fInc[2]==1) { if (fZb[fSlices[2]+1]>=xptnew) break; } else { if (fZb[fSlices[2]]<= xptnew) break; } } // printf(" %i/%i\n", fSlices[2], fIbz-2); } } // new indices are set -> Get new candidates if (fPriority[0]==1) { // we are entering the unique slice on this axis //---> intersect and store Y and Z if (fPriority[1]==2) { ndd[0] = fIndY[fOBy[fSlices[1]]]; if (!ndd[0]) return fCheckList; slice1 = (UChar_t*)(&fIndY[fOBy[fSlices[1]]+1]); islices++; } if (fPriority[2]==2) { ndd[1] = fIndZ[fOBz[fSlices[2]]]; if (!ndd[1]) return fCheckList; islices++; if (slice1) { slice2 = (UChar_t*)(&fIndZ[fOBz[fSlices[2]]+1]); } else { slice1 = (UChar_t*)(&fIndZ[fOBz[fSlices[2]]+1]); ndd[0] = ndd[1]; } } if (islices==1) { IntersectAndStore(ndd[0], slice1); } else { IntersectAndStore(ndd[0], slice1, ndd[1], slice2); } ncheck = fNcandidates; return fCheckList; } // We got into a new slice -> Get only new candidates left = (fInc[0]>0)?kTRUE:kFALSE; new_list = GetExtraX(fSlices[0], left, ncheck); // printf(" New list on X : %i new candidates\n", ncheck); if (!ncheck) return fCheckList; if (fPriority[1]==2) { ndd[0] = fIndY[fOBy[fSlices[1]]]; if (!ndd[0]) { ncheck = 0; return fCheckList; } slice1 = (UChar_t*)(&fIndY[fOBy[fSlices[1]]+1]); islices++; } if (fPriority[2]==2) { ndd[1] = fIndZ[fOBz[fSlices[2]]]; if (!ndd[1]) { ncheck = 0; return fCheckList; } islices++; if (slice1) { slice2 = (UChar_t*)(&fIndZ[fOBz[fSlices[2]]+1]); } else { slice1 = (UChar_t*)(&fIndZ[fOBz[fSlices[2]]+1]); ndd[0] = ndd[1]; } } if (!islices) return GetValidExtra(new_list, ncheck); if (islices==1) { return GetValidExtra(ndd[0], slice1, new_list, ncheck); } else { return GetValidExtra(ndd[0], slice1, ndd[1], slice2, new_list, ncheck); } case 1: if (isYlimit) return 0; // increment/decrement Y slice fSlices[1]=dind[1]; if (iforced) { // we have to recompute X and Z slices if (dslice>fLimits[0]) return 0; if (dslice>fLimits[2]) return 0; if ((dslice>dmin[0]) && fInc[0]) { xptnew = point[0]+dslice/fInvdir[0]; // printf(" recomputing X slice, pos=%g\n", xptnew); while (1) { fSlices[0] += fInc[0]; if (fInc[0]==1) { if (fXb[fSlices[0]+1]>=xptnew) break; } else { if (fXb[fSlices[0]]<= xptnew) break; } } // printf(" %i/%i\n", fSlices[0], fIbx-2); } if ((dslice>dmin[2]) && fInc[2]) { xptnew = point[2]+dslice/fInvdir[2]; // printf(" recomputing Z slice, pos=%g\n", xptnew); while (1) { fSlices[2] += fInc[2]; if (fInc[2]==1) { if (fZb[fSlices[2]+1]>=xptnew) break; } else { if (fZb[fSlices[2]]<= xptnew) break; } } // printf(" %i/%i\n", fSlices[2], fIbz-2); } } // new indices are set -> Get new candidates if (fPriority[1]==1) { // we are entering the unique slice on this axis //---> intersect and store X and Z if (fPriority[0]==2) { ndd[0] = fIndX[fOBx[fSlices[0]]]; if (!ndd[0]) return fCheckList; slice1 = (UChar_t*)(&fIndX[fOBx[fSlices[0]]+1]); islices++; } if (fPriority[2]==2) { ndd[1] = fIndZ[fOBz[fSlices[2]]]; if (!ndd[1]) return fCheckList; islices++; if (slice1) { slice2 = (UChar_t*)(&fIndZ[fOBz[fSlices[2]]+1]); } else { slice1 = (UChar_t*)(&fIndZ[fOBz[fSlices[2]]+1]); ndd[0] = ndd[1]; } } if (islices==1) { IntersectAndStore(ndd[0], slice1); } else { IntersectAndStore(ndd[0], slice1, ndd[1], slice2); } ncheck = fNcandidates; return fCheckList; } // We got into a new slice -> Get only new candidates left = (fInc[1]>0)?kTRUE:kFALSE; new_list = GetExtraY(fSlices[1], left, ncheck); // printf(" New list on Y : %i new candidates\n", ncheck); if (!ncheck) return fCheckList; if (fPriority[0]==2) { ndd[0] = fIndX[fOBx[fSlices[0]]]; if (!ndd[0]) { ncheck = 0; return fCheckList; } slice1 = (UChar_t*)(&fIndX[fOBx[fSlices[0]]+1]); islices++; } if (fPriority[2]==2) { ndd[1] = fIndZ[fOBz[fSlices[2]]]; if (!ndd[1]) { ncheck = 0; return fCheckList; } islices++; if (slice1) { slice2 = (UChar_t*)(&fIndZ[fOBz[fSlices[2]]+1]); } else { slice1 = (UChar_t*)(&fIndZ[fOBz[fSlices[2]]+1]); ndd[0] = ndd[1]; } } if (!islices) return GetValidExtra(new_list, ncheck); if (islices==1) { return GetValidExtra(ndd[0], slice1, new_list, ncheck); } else { return GetValidExtra(ndd[0], slice1, ndd[1], slice2, new_list, ncheck); } case 2: if (isZlimit) return 0; // increment/decrement Z slice fSlices[2]=dind[2]; if (iforced) { // we have to recompute Y and X slices if (dslice>fLimits[1]) return 0; if (dslice>fLimits[0]) return 0; if ((dslice>dmin[1]) && fInc[1]) { xptnew = point[1]+dslice/fInvdir[1]; // printf(" recomputing Y slice, pos=%g\n", xptnew); while (1) { fSlices[1] += fInc[1]; if (fInc[1]==1) { if (fYb[fSlices[1]+1]>=xptnew) break; } else { if (fYb[fSlices[1]]<= xptnew) break; } } // printf(" %i/%i\n", fSlices[1], fIby-2); } if ((dslice>dmin[0]) && fInc[0]) { xptnew = point[0]+dslice/fInvdir[0]; // printf(" recomputing X slice, pos=%g\n", xptnew); while (1) { fSlices[0] += fInc[0]; if (fInc[0]==1) { if (fXb[fSlices[0]+1]>=xptnew) break; } else { if (fXb[fSlices[0]]<= xptnew) break; } } // printf(" %i/%i\n", fSlices[0], fIbx-2); } } // new indices are set -> Get new candidates if (fPriority[2]==1) { // we are entering the unique slice on this axis //---> intersect and store Y and X if (fPriority[1]==2) { ndd[0] = fIndY[fOBy[fSlices[1]]]; if (!ndd[0]) return fCheckList; slice1 = (UChar_t*)(&fIndY[fOBy[fSlices[1]]+1]); islices++; } if (fPriority[0]==2) { ndd[1] = fIndX[fOBx[fSlices[0]]]; if (!ndd[1]) return fCheckList; islices++; if (slice1) { slice2 = (UChar_t*)(&fIndX[fOBx[fSlices[0]]+1]); } else { slice1 = (UChar_t*)(&fIndX[fOBx[fSlices[0]]+1]); ndd[0] = ndd[1]; } } if (islices==1) { IntersectAndStore(ndd[0], slice1); } else { IntersectAndStore(ndd[0], slice1, ndd[1], slice2); } ncheck = fNcandidates; return fCheckList; } // We got into a new slice -> Get only new candidates left = (fInc[2]>0)?kTRUE:kFALSE; new_list = GetExtraZ(fSlices[2], left, ncheck); // printf(" New list on Z : %i new candidates\n", ncheck); if (!ncheck) return fCheckList; if (fPriority[1]==2) { ndd[0] = fIndY[fOBy[fSlices[1]]]; if (!ndd[0]) { ncheck = 0; return fCheckList; } slice1 = (UChar_t*)(&fIndY[fOBy[fSlices[1]]+1]); islices++; } if (fPriority[0]==2) { ndd[1] = fIndX[fOBx[fSlices[0]]]; if (!ndd[1]) { ncheck = 0; return fCheckList; } islices++; if (slice1) { slice2 = (UChar_t*)(&fIndX[fOBx[fSlices[0]]+1]); } else { slice1 = (UChar_t*)(&fIndX[fOBx[fSlices[0]]+1]); ndd[0] = ndd[1]; } } if (!islices) return GetValidExtra(new_list, ncheck); if (islices==1) { return GetValidExtra(ndd[0], slice1, new_list, ncheck); } else { return GetValidExtra(ndd[0], slice1, ndd[1], slice2, new_list, ncheck); } default: Error("GetNextCandidates", "Invalid islice=%i inside %s", islice, fVolume->GetName()); } return 0; } //----------------------------------------------------------------------------- void TGeoVoxelFinder::SortCrossedVoxels(Double_t *point, Double_t *dir) { // get the list in the next voxel crossed by a ray fCurrentVoxel = 0; // printf("###Sort crossed voxels for %s\n", fVolume->GetName()); fNcandidates = 0; Int_t loc = 1+((fVolume->GetNdaughters()-1)>>3); // printf(" LOC=%i\n", loc*sizeof(UChar_t)); // UChar_t *bits = gGeoManager->GetBits(); memset(fBits1, 0, loc); memset(fInc, 0, 3*sizeof(Int_t)); for (Int_t i=0; i<3; i++) { fInvdir[i] = TGeoShape::Big(); if (TMath::Abs(dir[i])<1E-10) continue; fInc[i] = (dir[i]>0)?1:-1; fInvdir[i] = 1./dir[i]; } Bool_t flag = GetIndices(point); TGeoBBox *box = (TGeoBBox*)(fVolume->GetShape()); if (fInc[0]==0) { fLimits[0] = TGeoShape::Big(); } else { if (fSlices[0]==-2) { // no slice on this axis -> get limit to bounding box limit fLimits[0] = ((box->GetOrigin())[0]-point[0]+fInc[0]*box->GetDX())*fInvdir[0]; } else { if (fInc[0]==1) { fLimits[0] = (fXb[fIbx-1]-point[0])*fInvdir[0]; } else { fLimits[0] = (fXb[0]-point[0])*fInvdir[0]; } } } if (fInc[1]==0) { fLimits[1] = TGeoShape::Big(); } else { if (fSlices[1]==-2) { // no slice on this axis -> get limit to bounding box limit fLimits[1] = ((box->GetOrigin())[1]-point[1]+fInc[1]*box->GetDY())*fInvdir[1]; } else { if (fInc[1]==1) { fLimits[1] = (fYb[fIby-1]-point[1])*fInvdir[1]; } else { fLimits[1] = (fYb[0]-point[1])*fInvdir[1]; } } } if (fInc[2]==0) { fLimits[2] = TGeoShape::Big(); } else { if (fSlices[2]==-2) { // no slice on this axis -> get limit to bounding box limit fLimits[2] = ((box->GetOrigin())[2]-point[2]+fInc[2]*box->GetDZ())*fInvdir[2]; } else { if (fInc[2]==1) { fLimits[2] = (fZb[fIbz-1]-point[2])*fInvdir[2]; } else { fLimits[2] = (fZb[0]-point[2])*fInvdir[2]; } } } if (!flag) { // printf(" NO candidates in first voxel\n"); // printf(" bits[0]=%i\n", bits[0]); return; } // printf(" current slices : %i %i %i\n", fSlices[0], fSlices[1], fSlices[2]); Int_t nd[3]; Int_t islices = 0; memset(&nd[0], 0, 3*sizeof(Int_t)); UChar_t *slicex = 0; if (fPriority[0]==2) { nd[0] = fIndX[fOBx[fSlices[0]]]; slicex=(UChar_t*)(&fIndX[fOBx[fSlices[0]]+1]); islices++; } UChar_t *slicey = 0; if (fPriority[1]==2) { nd[1] = fIndY[fOBy[fSlices[1]]]; islices++; if (slicex) { slicey=(UChar_t*)(&fIndY[fOBy[fSlices[1]]+1]); } else { slicex=(UChar_t*)(&fIndY[fOBy[fSlices[1]]+1]); nd[0] = nd[1]; } } UChar_t *slicez = 0; if (fPriority[2]==2) { nd[2] = fIndZ[fOBz[fSlices[2]]]; islices++; if (slicex && slicey) { slicez=(UChar_t*)(&fIndZ[fOBz[fSlices[2]]+1]); } else { if (slicex) { slicey=(UChar_t*)(&fIndZ[fOBz[fSlices[2]]+1]); nd[1] = nd[2]; } else { slicex=(UChar_t*)(&fIndZ[fOBz[fSlices[2]]+1]); nd[0] = nd[2]; } } } // printf("Ndaughters in first voxel : %i %i %i\n", nd[0], nd[1], nd[2]); switch (islices) { case 0: Error("SortCrossedVoxels", "no slices for %s", fVolume->GetName()); // printf("Slices :(%i,%i,%i) Priority:(%i,%i,%i)\n", fSlices[0], fSlices[1], fSlices[2], fPriority[0], fPriority[1], fPriority[2]); return; case 1: IntersectAndStore(nd[0], slicex); break; case 2: IntersectAndStore(nd[0], slicex, nd[1], slicey); break; default: IntersectAndStore(nd[0], slicex, nd[1], slicey, nd[2], slicez); } // printf(" bits[0]=%i END\n", bits[0]); // if (fNcandidates) { // printf(" candidates for first voxel :\n"); // for (Int_t i=0; i<fNcandidates; i++) printf(" %i\n", fCheckList[i]); // } } //----------------------------------------------------------------------------- Int_t *TGeoVoxelFinder::GetCheckList(Double_t *point, Int_t &nelem) { // get the list of daughter indices for which point is inside their bbox // if (!fBoxes) return 0; if (fVolume->GetNdaughters() == 1) { if (fXb) { if (point[0]<fXb[0] || point[0]>fXb[1]) return 0; } if (fYb) { if (point[1]<fYb[0] || point[1]>fYb[1]) return 0; } if (fZb) { if (point[2]<fZb[0] || point[2]>fZb[1]) return 0; } fCheckList[0] = 0; nelem = 1; return fCheckList; } Int_t nslices = 0; UChar_t *slice1 = 0; UChar_t *slice2 = 0; UChar_t *slice3 = 0; Int_t nd[3] = {0,0,0}; Int_t im; if (fPriority[0]) { im = TMath::BinarySearch(fIbx, fXb, point[0]); if ((im==-1) || (im==fIbx-1)) return 0; if (fPriority[0]==2) { nd[0] = fIndX[fOBx[im]]; if (!nd[0]) return 0; nslices++; slice1 = (UChar_t*)(&fIndX[fOBx[im]+1]); } } if (fPriority[1]) { im = TMath::BinarySearch(fIby, fYb, point[1]); if ((im==-1) || (im==fIby-1)) return 0; if (fPriority[1]==2) { nd[1] = fIndY[fOBy[im]]; if (!nd[1]) return 0; nslices++; if (slice1) { slice2 = (UChar_t*)(&fIndY[fOBy[im]+1]); } else { slice1 = (UChar_t*)(&fIndY[fOBy[im]+1]); nd[0] = nd[1]; } } } if (fPriority[2]) { im = TMath::BinarySearch(fIbz, fZb, point[2]); if ((im==-1) || (im==fIbz-1)) return 0; if (fPriority[2]==2) { nd[2] = fIndZ[fOBz[im]]; if (!nd[2]) return 0; nslices++; if (slice1 && slice2) { slice3 = (UChar_t*)(&fIndZ[fOBz[im]+1]); } else { if (slice1) { slice2 = (UChar_t*)(&fIndZ[fOBz[im]+1]); nd[1] = nd[2]; } else { slice1 = (UChar_t*)(&fIndZ[fOBz[im]+1]); nd[0] = nd[2]; } } } } nelem = 0; // Int_t i = 0; Bool_t intersect = kFALSE; switch (nslices) { case 0: Error("GetCheckList", "No slices for %s", fVolume->GetName()); return 0; case 1: intersect = Intersect(nd[0], slice1, nelem, fCheckList); break; case 2: intersect = Intersect(nd[0], slice1, nd[1], slice2, nelem, fCheckList); break; default: intersect = Intersect(nd[0], slice1, nd[1], slice2, nd[2], slice3, nelem, fCheckList); } if (intersect) return fCheckList; return 0; } //----------------------------------------------------------------------------- Int_t *TGeoVoxelFinder::GetVoxelCandidates(Int_t i, Int_t j, Int_t k, Int_t &ncheck) { // get the list of candidates in voxel (i,j,k) - no check UChar_t *slice1 = 0; UChar_t *slice2 = 0; UChar_t *slice3 = 0; Int_t nd[3] = {0,0,0}; Int_t nslices = 0; if (fPriority[0]==2) { nd[0] = fIndX[fOBx[i]]; if (!nd[0]) return 0; nslices++; slice1 = (UChar_t*)(&fIndX[fOBx[i]+1]); } if (fPriority[1]==2) { nd[1] = fIndY[fOBy[j]]; if (!nd[1]) return 0; nslices++; if (slice1) { slice2 = (UChar_t*)(&fIndY[fOBy[j]+1]); } else { slice1 = (UChar_t*)(&fIndY[fOBy[j]+1]); nd[0] = nd[1]; } } if (fPriority[2]==2) { nd[2] = fIndZ[fOBz[k]]; if (!nd[2]) return 0; nslices++; if (slice1 && slice2) { slice3 = (UChar_t*)(&fIndZ[fOBz[k]+1]); } else { if (slice1) { slice2 = (UChar_t*)(&fIndZ[fOBz[k]+1]); nd[1] = nd[2]; } else { slice1 = (UChar_t*)(&fIndZ[fOBz[k]+1]); nd[0] = nd[2]; } } } Bool_t intersect = kFALSE; switch (nslices) { case 0: Error("GetCheckList", "No slices for %s", fVolume->GetName()); return 0; case 1: intersect = Intersect(nd[0], slice1, ncheck, fCheckList); break; case 2: intersect = Intersect(nd[0], slice1, nd[1], slice2, ncheck, fCheckList); break; default: intersect = Intersect(nd[0], slice1, nd[1], slice2, nd[2], slice3, ncheck, fCheckList); } if (intersect) return fCheckList; return 0; } //----------------------------------------------------------------------------- Int_t *TGeoVoxelFinder::GetNextVoxel(Double_t *point, Double_t * /*dir*/, Int_t &ncheck) { // get the list of new candidates for the next voxel crossed by current ray // printf("### GetNextVoxel\n"); if (fCurrentVoxel==0) { // printf(">>> first voxel, %i candidates\n", ncheck); // printf(" bits[0]=%i\n", gGeoManager->GetBits()[0]); fCurrentVoxel++; ncheck = fNcandidates; return fCheckList; } fCurrentVoxel++; // printf(">>> voxel %i\n", fCurrentVoxel); // Get slices for next voxel // printf("before - fSlices : %i %i %i\n", fSlices[0], fSlices[1], fSlices[2]); return GetNextCandidates(point, ncheck); } //----------------------------------------------------------------------------- Bool_t TGeoVoxelFinder::Intersect(Int_t n1, UChar_t *array1, Int_t &nf, Int_t *result) { // return the list of nodes corresponding to one array of bits Int_t nd = fVolume->GetNdaughters(); // also number of bits to scan nf = 0; Int_t nbytes = 1+((nd-1)>>3); Int_t current_byte; Int_t current_bit; UChar_t byte; Bool_t ibreak = kFALSE; for (current_byte=0; current_byte<nbytes; current_byte++) { byte = array1[current_byte]; if (!byte) continue; for (current_bit=0; current_bit<8; current_bit++) { if (byte & (1<<current_bit)) { result[nf++] = (current_byte<<3)+current_bit; if (nf==n1) { ibreak = kTRUE; break; } } } if (ibreak) return kTRUE; } return kTRUE; } //----------------------------------------------------------------------------- Bool_t TGeoVoxelFinder::IntersectAndStore(Int_t n1, UChar_t *array1) { // return the list of nodes corresponding to one array of bits Int_t nd = fVolume->GetNdaughters(); // also number of bits to scan // UChar_t *bits = gGeoManager->GetBits(); fNcandidates = 0; Int_t nbytes = 1+((nd-1)>>3); memcpy(fBits1, array1, nbytes*sizeof(UChar_t)); Int_t current_byte; Int_t current_bit; UChar_t byte; Bool_t ibreak = kFALSE; Int_t icand; for (current_byte=0; current_byte<nbytes; current_byte++) { byte = array1[current_byte]; icand = current_byte<<3; if (!byte) continue; for (current_bit=0; current_bit<8; current_bit++) { if (byte & (1<<current_bit)) { fCheckList[fNcandidates++] = icand+current_bit; if (fNcandidates==n1) { ibreak = kTRUE; break; } } } if (ibreak) return kTRUE; } return kTRUE; } //----------------------------------------------------------------------------- Bool_t TGeoVoxelFinder::Union(Int_t n1, UChar_t *array1) { // make union of older bits with new array // printf("Union - one slice\n"); Int_t nd = fVolume->GetNdaughters(); // also number of bits to scan // UChar_t *bits = gGeoManager->GetBits(); fNcandidates = 0; Int_t nbytes = 1+((nd-1)>>3); Int_t current_byte; Int_t current_bit; UChar_t byte; Bool_t ibreak = kFALSE; for (current_byte=0; current_byte<nbytes; current_byte++) { // printf(" byte %i : bits=%i array=%i\n", current_byte, bits[current_byte], array1[current_byte]); byte = (~fBits1[current_byte]) & array1[current_byte]; if (!byte) continue; for (current_bit=0; current_bit<8; current_bit++) { if (byte & (1<<current_bit)) { fCheckList[fNcandidates++] = (current_byte<<3)+current_bit; if (fNcandidates==n1) { ibreak = kTRUE; break; } } } fBits1[current_byte] |= byte; if (ibreak) return kTRUE; } return (fNcandidates>0); } //----------------------------------------------------------------------------- Bool_t TGeoVoxelFinder::Union(Int_t /*n1*/, UChar_t *array1, Int_t /*n2*/, UChar_t *array2) { // make union of older bits with new array // printf("Union - two slices\n"); Int_t nd = fVolume->GetNdaughters(); // also number of bits to scan // UChar_t *bits = gGeoManager->GetBits(); fNcandidates = 0; Int_t nbytes = 1+((nd-1)>>3); Int_t current_byte; Int_t current_bit; UChar_t byte; for (current_byte=0; current_byte<nbytes; current_byte++) { byte = (~fBits1[current_byte]) & (array1[current_byte] & array2[current_byte]); if (!byte) continue; for (current_bit=0; current_bit<8; current_bit++) { if (byte & (1<<current_bit)) { fCheckList[fNcandidates++] = (current_byte<<3)+current_bit; } } fBits1[current_byte] |= byte; } return (fNcandidates>0); } //----------------------------------------------------------------------------- Bool_t TGeoVoxelFinder::Union(Int_t /*n1*/, UChar_t *array1, Int_t /*n2*/, UChar_t *array2, Int_t /*n3*/, UChar_t *array3) { // make union of older bits with new array // printf("Union - three slices\n"); // printf("n1=%i n2=%i n3=%i\n", n1,n2,n3); Int_t nd = fVolume->GetNdaughters(); // also number of bits to scan // UChar_t *bits = gGeoManager->GetBits(); fNcandidates = 0; Int_t nbytes = 1+((nd-1)>>3); Int_t current_byte; Int_t current_bit; UChar_t byte; for (current_byte=0; current_byte<nbytes; current_byte++) { byte = (~fBits1[current_byte]) & (array1[current_byte] & array2[current_byte] & array3[current_byte]); if (!byte) continue; for (current_bit=0; current_bit<8; current_bit++) { if (byte & (1<<current_bit)) { fCheckList[fNcandidates++] = (current_byte<<3)+current_bit; } } fBits1[current_byte] |= byte; } return (fNcandidates>0); } //----------------------------------------------------------------------------- Bool_t TGeoVoxelFinder::Intersect(Int_t n1, UChar_t *array1, Int_t n2, UChar_t *array2, Int_t &nf, Int_t *result) { // return the list of nodes corresponding to the intersection of two arrays of bits Int_t nd = fVolume->GetNdaughters(); // also number of bits to scan nf = 0; Int_t nbytes = 1+((nd-1)>>3); Int_t current_byte; Int_t current_bit; UChar_t byte; Bool_t ibreak = kFALSE; for (current_byte=0; current_byte<nbytes; current_byte++) { byte = array1[current_byte] & array2[current_byte]; if (!byte) continue; for (current_bit=0; current_bit<8; current_bit++) { if (byte & (1<<current_bit)) { result[nf++] = (current_byte<<3)+current_bit; if ((nf==n1) || (nf==n2)) { ibreak = kTRUE; break; } } } if (ibreak) return kTRUE; } return (nf>0); } //----------------------------------------------------------------------------- Bool_t TGeoVoxelFinder::IntersectAndStore(Int_t /*n1*/, UChar_t *array1, Int_t /*n2*/, UChar_t *array2) { // return the list of nodes corresponding to the intersection of two arrays of bits Int_t nd = fVolume->GetNdaughters(); // also number of bits to scan // UChar_t *bits = gGeoManager->GetBits(); fNcandidates = 0; Int_t nbytes = 1+((nd-1)>>3); // memset(bits, 0, nbytes*sizeof(UChar_t)); Int_t current_byte; Int_t current_bit; Int_t icand; UChar_t byte; for (current_byte=0; current_byte<nbytes; current_byte++) { byte = array1[current_byte] & array2[current_byte]; icand = current_byte<<3; fBits1[current_byte] = byte; if (!byte) continue; for (current_bit=0; current_bit<8; current_bit++) { if (byte & (1<<current_bit)) { fCheckList[fNcandidates++] = icand+current_bit; } } } return (fNcandidates>0); } //----------------------------------------------------------------------------- Bool_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) { // return the list of nodes corresponding to the intersection of three arrays of bits Int_t nd = fVolume->GetNdaughters(); // also number of bits to scan nf = 0; Int_t nbytes = 1+((nd-1)>>3); Int_t current_byte; Int_t current_bit; UChar_t byte; Bool_t ibreak = kFALSE; for (current_byte=0; current_byte<nbytes; current_byte++) { byte = array1[current_byte] & array2[current_byte] & array3[current_byte]; if (!byte) continue; for (current_bit=0; current_bit<8; current_bit++) { if (byte & (1<<current_bit)) { result[nf++] = (current_byte<<3)+current_bit; if ((nf==n1) || (nf==n2) || (nf==n3)) { ibreak = kTRUE; break; } } } if (ibreak) return kTRUE; } return (nf>0); } //----------------------------------------------------------------------------- Bool_t TGeoVoxelFinder::IntersectAndStore(Int_t /*n1*/, UChar_t *array1, Int_t /*n2*/, UChar_t *array2, Int_t /*n3*/, UChar_t *array3) { // return the list of nodes corresponding to the intersection of three arrays of bits Int_t nd = fVolume->GetNdaughters(); // also number of bits to scan // UChar_t *bits = gGeoManager->GetBits(); fNcandidates = 0; Int_t nbytes = 1+((nd-1)>>3); // memset(bits, 0, nbytes*sizeof(UChar_t)); Int_t current_byte; Int_t current_bit; Int_t icand; UChar_t byte; for (current_byte=0; current_byte<nbytes; current_byte++) { byte = array1[current_byte] & array2[current_byte] & array3[current_byte]; icand = current_byte<<3; fBits1[current_byte] = byte; if (!byte) continue; for (current_bit=0; current_bit<8; current_bit++) { if (byte & (1<<current_bit)) { fCheckList[fNcandidates++] = icand+current_bit; } } } return (fNcandidates>0); } //----------------------------------------------------------------------------- void TGeoVoxelFinder::SortAll(Option_t *) { // order bounding boxes along x, y, z Int_t nd = fVolume->GetNdaughters(); Int_t nperslice = 1 /*N in slice*/ + 1+(nd-1)/(8*sizeof(Int_t)); /*Nbytes per slice*/ Int_t nmaxslices = 2*nd-1; // max number of slices on each axis Double_t *boundaries = new Double_t[6*nd]; // list of different boundaries Double_t xmin, xmax, ymin, ymax, zmin, zmax; TGeoBBox *box = (TGeoBBox*)fVolume->GetShape(); // bounding box for volume // compute range on X, Y, Z according to volume bounding box xmin = (box->GetOrigin())[0] - box->GetDX(); xmax = (box->GetOrigin())[0] + box->GetDX(); ymin = (box->GetOrigin())[1] - box->GetDY(); ymax = (box->GetOrigin())[1] + box->GetDY(); zmin = (box->GetOrigin())[2] - box->GetDZ(); zmax = (box->GetOrigin())[2] + box->GetDZ(); if ((xmin>=xmax) || (ymin>=ymax) || (zmin>=zmax)) { Error("SortAll", "Wrong bounding box for volume %s", fVolume->GetName()); return; } Int_t id; // compute boundaries coordinates on X,Y,Z for (id=0; id<nd; id++) { // x boundaries boundaries[2*id] = fBoxes[6*id+3]-fBoxes[6*id]; boundaries[2*id+1] = fBoxes[6*id+3]+fBoxes[6*id]; // y boundaries boundaries[2*id+2*nd] = fBoxes[6*id+4]-fBoxes[6*id+1]; boundaries[2*id+2*nd+1] = fBoxes[6*id+4]+fBoxes[6*id+1]; // z boundaries boundaries[2*id+4*nd] = fBoxes[6*id+5]-fBoxes[6*id+2]; boundaries[2*id+4*nd+1] = fBoxes[6*id+5]+fBoxes[6*id+2]; } Int_t *index = new Int_t[2*nd]; // indexes for sorted boundaries on one axis Int_t *ind = new Int_t[nmaxslices*nperslice]; // ind[fOBx[i]] = ndghts in slice fInd[i]--fInd[i+1] Int_t *extra = new Int_t[nmaxslices*4]; // extra[fOEx[i]] = nextra_to_left (i/i-1) // extra[fOEx[i]+1] = nextra_to_right (i/i+1) // Int_t *extra_to_left = extra[fOEx[i]+2] // Int_t *extra_to_right = extra[fOEx[i]+2+nextra_to_left] Double_t *temp = new Double_t[2*nd]; // temporary array to store sorted boundary positions Int_t current = 0; Int_t indextra = 0; Int_t nleft, nright; Int_t *extra_left = new Int_t[nd]; Int_t *extra_right = new Int_t[nd]; Double_t xxmin, xxmax, xbmin, xbmax, ddx1, ddx2; UChar_t *bits; UInt_t loc, bitnumber; UChar_t bit; // sort x boundaries Int_t ib = 0; // total number of DIFFERENT boundaries TMath::Sort(2*nd, &boundaries[0], &index[0], kFALSE); // compact common boundaries for (id=0; id<2*nd; id++) { if (!ib) {temp[ib++] = boundaries[index[id]]; continue;} if (TMath::Abs(temp[ib-1]-boundaries[index[id]])>1E-10) temp[ib++] = boundaries[index[id]]; } // now find priority if (ib < 2) { Error("SortAll", "Cannot voxelize %s :less than 2 boundaries on X", fVolume->GetName()); delete [] index; delete [] ind; delete [] temp; delete [] extra; delete [] extra_left; delete [] extra_right; SetInvalid(); return; } if (ib == 2) { // check range if (((temp[0]-xmin)<1E-10) && ((temp[1]-xmax)>-1E-10)) { // ordering on this axis makes no sense. Clear all arrays. fPriority[0] = 0; // always skip this axis if (fIndX) delete [] fIndX; fIndX = 0; fNx = 0; if (fXb) delete [] fXb; fXb = 0; fIbx = 0; if (fOBx) delete [] fOBx; fOBx = 0; fNox = 0; } else { fPriority[0] = 1; // all in one slice } } else { fPriority[0] = 2; // check all } // store compacted boundaries if (fPriority[0]) { if (fXb) delete [] fXb; fXb = new Double_t[ib]; memcpy(fXb, &temp[0], ib*sizeof(Double_t)); fIbx = ib; } //--- now build the lists of nodes in each slice of X axis if (fPriority[0]==2) { memset(ind, 0, (nmaxslices*nperslice)*sizeof(Int_t)); if (fOBx) delete [] fOBx; fNox = fIbx-1; // number of different slices fOBx = new Int_t[fNox]; // offsets in ind if (fOEx) delete [] fOEx; fOEx = new Int_t[fNox]; // offsets in extra current = 0; indextra = 0; //--- now loop all slices for (id=0; id<fNox; id++) { fOBx[id] = current; // offset in dght list for this slice fOEx[id] = indextra; // offset in exta list for this slice ind[current] = 0; // ndght in this slice extra[indextra] = extra[indextra+1] = 0; // no extra left/right nleft = nright = 0; bits = (UChar_t*)(&ind[current+1]); // adress of bits for this slice xxmin = fXb[id]; xxmax = fXb[id+1]; for (Int_t ic=0; ic<nd; ic++) { xbmin = fBoxes[6*ic+3]-fBoxes[6*ic]; xbmax = fBoxes[6*ic+3]+fBoxes[6*ic]; ddx1 = xbmin-xxmax; if (ddx1>-1E-10) continue; ddx2 = xbmax-xxmin; if (ddx2<1E-10) continue; // daughter ic in interval //---> set the bit ind[current]++; bitnumber = (UInt_t)ic; loc = bitnumber/8; bit = bitnumber%8; bits[loc] |= 1<<bit; //---> chech if it is extra to left/right //--- left ddx1 = xbmin-xxmin; ddx2 = xbmax-xxmax; if ((id==0) || (ddx1>-1E-8)) { extra_left[nleft++] = ic; } //---right if ((id==(fNoz-1)) || (ddx2<1E-8)) { extra_right[nright++] = ic; } } //--- compute offset of next slice if (ind[current]>0) { current += nperslice; } else { current ++; } //--- copy extra candidates extra[indextra] = nleft; extra[indextra+1] = nright; if (nleft) memcpy(&extra[indextra+2], extra_left, nleft*sizeof(Int_t)); if (nright) memcpy(&extra[indextra+2+nleft], extra_right, nright*sizeof(Int_t)); indextra += 2+nleft+nright; } if (fIndX) delete [] fIndX; fNx = current; fIndX = new Int_t[current]; memcpy(fIndX, ind, current*sizeof(Int_t)); if (fExtraX) delete [] fExtraX; fNex = indextra; if (indextra>nmaxslices*4) printf("Woops!!!\n"); fExtraX = new Int_t[indextra]; memcpy(fExtraX, extra, indextra*sizeof(Int_t)); } // sort y boundaries ib = 0; TMath::Sort(2*nd, &boundaries[2*nd], &index[0], kFALSE); // compact common boundaries for (id=0; id<2*nd; id++) { if (!ib) {temp[ib++] = boundaries[2*nd+index[id]]; continue;} if (TMath::Abs(temp[ib-1]-boundaries[2*nd+index[id]])>1E-10) temp[ib++]=boundaries[2*nd+index[id]]; } // now find priority on Y if (ib < 2) { Error("SortAll", "Cannot voxelize %s :less than 2 boundaries on Y", fVolume->GetName()); delete [] index; delete [] ind; delete [] temp; delete [] extra; delete [] extra_left; delete [] extra_right; SetInvalid(); return; } if (ib == 2) { // check range if (((temp[0]-ymin)<1E-10) && ((temp[1]-ymax)>-1E-10)) { // ordering on this axis makes no sense. Clear all arrays. fPriority[1] = 0; // always skip this axis if (fIndY) delete [] fIndY; fIndY = 0; fNy = 0; if (fYb) delete [] fYb; fYb = 0; fIby = 0; if (fOBy) delete [] fOBy; fOBy = 0; fNoy = 0; } else { fPriority[1] = 1; // all in one slice } } else { fPriority[1] = 2; // check all } // store compacted boundaries if (fPriority[1]) { if (fYb) delete [] fYb; fYb = new Double_t[ib]; memcpy(fYb, &temp[0], ib*sizeof(Double_t)); fIby = ib; } if (fPriority[1]==2) { //--- now build the lists of nodes in each slice of Y axis memset(ind, 0, (nmaxslices*nperslice)*sizeof(Int_t)); if (fOBy) delete [] fOBy; fNoy = fIby-1; // number of different slices fOBy = new Int_t[fNoy]; // offsets in ind if (fOEy) delete [] fOEy; fOEy = new Int_t[fNoy]; // offsets in extra current = 0; indextra = 0; //--- now loop all slices for (id=0; id<fNoy; id++) { fOBy[id] = current; // offset of dght list fOEy[id] = indextra; // offset in exta list for this slice ind[current] = 0; // ndght in this slice extra[indextra] = extra[indextra+1] = 0; // no extra left/right nleft = nright = 0; bits = (UChar_t*)(&ind[current+1]); // adress of bits for this slice xxmin = fYb[id]; xxmax = fYb[id+1]; for (Int_t ic=0; ic<nd; ic++) { xbmin = fBoxes[6*ic+4]-fBoxes[6*ic+1]; xbmax = fBoxes[6*ic+4]+fBoxes[6*ic+1]; ddx1 = xbmin-xxmax; if (ddx1>-1E-10) continue; ddx2 = xbmax-xxmin; if (ddx2<1E-10) continue; // daughter ic in interval //---> set the bit ind[current]++; bitnumber = (UInt_t)ic; loc = bitnumber/8; bit = bitnumber%8; bits[loc] |= 1<<bit; //---> chech if it is extra to left/right //--- left ddx1 = xbmin-xxmin; ddx2 = xbmax-xxmax; if ((id==0) || (ddx1>-1E-8)) { extra_left[nleft++] = ic; } //---right if ((id==(fNoz-1)) || (ddx2<1E-8)) { extra_right[nright++] = ic; } } //--- compute offset of next slice if (ind[current]>0) { current += nperslice; } else { current ++; } //--- copy extra candidates extra[indextra] = nleft; extra[indextra+1] = nright; if (nleft) memcpy(&extra[indextra+2], extra_left, nleft*sizeof(Int_t)); if (nright) memcpy(&extra[indextra+2+nleft], extra_right, nright*sizeof(Int_t)); indextra += 2+nleft+nright; } if (fIndY) delete [] fIndY; fNy = current; fIndY = new Int_t[current]; memcpy(fIndY, &ind[0], current*sizeof(Int_t)); if (fExtraY) delete [] fExtraY; fNey = indextra; if (indextra>nmaxslices*4) printf("Woops!!!\n"); fExtraY = new Int_t[indextra]; memcpy(fExtraY, extra, indextra*sizeof(Int_t)); } // sort z boundaries ib = 0; TMath::Sort(2*nd, &boundaries[4*nd], &index[0], kFALSE); // compact common boundaries for (id=0; id<2*nd; id++) { if (!ib) {temp[ib++] = boundaries[4*nd+index[id]]; continue;} if ((TMath::Abs(temp[ib-1]-boundaries[4*nd+index[id]]))>1E-10) temp[ib++]=boundaries[4*nd+index[id]]; } // now find priority on Z if (ib < 2) { Error("SortAll", "Cannot voxelize %s :less than 2 boundaries on Z", fVolume->GetName()); delete [] index; delete [] ind; delete [] temp; delete [] extra; delete [] extra_left; delete [] extra_right; SetInvalid(); return; } if (ib == 2) { // check range if (((temp[0]-zmin)<1E-10) && ((temp[1]-zmax)>-1E-10)) { // ordering on this axis makes no sense. Clear all arrays. fPriority[2] = 0; if (fIndZ) delete [] fIndZ; fIndZ = 0; fNz = 0; if (fZb) delete [] fZb; fZb = 0; fIbz = 0; if (fOBz) delete [] fOBz; fOBz = 0; fNoz = 0; } else { fPriority[2] = 1; // all in one slice } } else { fPriority[2] = 2; // check all } // store compacted boundaries if (fPriority[2]) { if (fZb) delete [] fZb; fZb = new Double_t[ib]; memcpy(fZb, &temp[0], ib*sizeof(Double_t)); fIbz = ib; } if (fPriority[2]==2) { //--- now build the lists of nodes in each slice of Y axis memset(ind, 0, (nmaxslices*nperslice)*sizeof(Int_t)); if (fOBz) delete [] fOBz; fNoz = fIbz-1; // number of different slices fOBz = new Int_t[fNoz]; // offsets in ind if (fOEz) delete [] fOEz; fOEz = new Int_t[fNoz]; // offsets in extra current = 0; indextra = 0; //--- now loop all slices for (id=0; id<fNoz; id++) { fOBz[id] = current; // offset of dght list fOEz[id] = indextra; // offset in exta list for this slice ind[current] = 0; // ndght in this slice extra[indextra] = extra[indextra+1] = 0; // no extra left/right nleft = nright = 0; bits = (UChar_t*)(&ind[current+1]); // adress of bits for this slice xxmin = fZb[id]; xxmax = fZb[id+1]; for (Int_t ic=0; ic<nd; ic++) { xbmin = fBoxes[6*ic+5]-fBoxes[6*ic+2]; xbmax = fBoxes[6*ic+5]+fBoxes[6*ic+2]; ddx1 = xbmin-xxmax; if (ddx1>-1E-10) continue; ddx2 = xbmax-xxmin; if (ddx2<1E-10) continue; // daughter ic in interval //---> set the bit ind[current]++; bitnumber = (UInt_t)ic; loc = bitnumber/8; bit = bitnumber%8; bits[loc] |= 1<<bit; //---> chech if it is extra to left/right //--- left ddx1 = xbmin-xxmin; ddx2 = xbmax-xxmax; if ((id==0) || (ddx1>-1E-8)) { extra_left[nleft++] = ic; } //---right if ((id==(fNoz-1)) || (ddx2<1E-8)) { extra_right[nright++] = ic; } } //--- compute offset of next slice if (ind[current]>0) { current += nperslice; } else { current ++; } //--- copy extra candidates extra[indextra] = nleft; extra[indextra+1] = nright; if (nleft) memcpy(&extra[indextra+2], extra_left, nleft*sizeof(Int_t)); if (nright) memcpy(&extra[indextra+2+nleft], extra_right, nright*sizeof(Int_t)); indextra += 2+nleft+nright; } if (fIndZ) delete [] fIndZ; fNz = current; fIndZ = new Int_t[current]; memcpy(fIndZ, &ind[0], current*sizeof(Int_t)); if (fExtraZ) delete [] fExtraZ; fNez = indextra; if (indextra>nmaxslices*4) printf("Woops!!!\n"); fExtraZ = new Int_t[indextra]; memcpy(fExtraZ, extra, indextra*sizeof(Int_t)); } delete [] boundaries; delete [] index; delete [] temp; delete [] ind; delete [] extra; delete [] extra_left; delete [] extra_right; // Print(); if ((!fPriority[0]) && (!fPriority[1]) && (!fPriority[2])) { fVolume->SetVoxelFinder(0); delete this; } else { // Efficiency(); } } //----------------------------------------------------------------------------- void TGeoVoxelFinder::Print(Option_t *) const { Int_t id, i; printf("Voxels for volume %s (nd=%i)\n", fVolume->GetName(), fVolume->GetNdaughters()); printf("priority : x=%i y=%i z=%i\n", fPriority[0], fPriority[1], fPriority[2]); // return; Int_t nextra; Int_t nbytes = 1+((fVolume->GetNdaughters()-1)>>3); UChar_t byte, bit; UChar_t *slice; printf("XXX\n"); if (fPriority[0]==2) { for (id=0; id<fIbx; id++) { printf("%15.10f\n",fXb[id]); if (id == (fIbx-1)) break; printf("slice %i : %i\n", id, fIndX[fOBx[id]]); if (fIndX[fOBx[id]]) { slice = (UChar_t*)(&fIndX[fOBx[id]+1]); for (i=0; i<nbytes; i++) { byte = slice[i]; for (bit=0; bit<8; bit++) { if (byte & (1<<bit)) printf(" %i ", 8*i+bit); } } printf("\n"); } GetExtraX(id,kTRUE,nextra); printf(" extra_about_left = %i\n", nextra); GetExtraX(id,kFALSE,nextra); printf(" extra_about_right = %i\n", nextra); } } printf("YYY\n"); if (fPriority[1]==2) { for (id=0; id<fIby; id++) { printf("%15.10f\n", fYb[id]); if (id == (fIby-1)) break; printf("slice %i : %i\n", id, fIndY[fOBy[id]]); if (fIndY[fOBy[id]]) { slice = (UChar_t*)(&fIndY[fOBy[id]+1]); for (i=0; i<nbytes; i++) { byte = slice[i]; for (bit=0; bit<8; bit++) { if (byte & (1<<bit)) printf(" %i ", 8*i+bit); } } } GetExtraY(id,kTRUE,nextra); printf(" extra_about_left = %i\n", nextra); GetExtraY(id,kFALSE,nextra); printf(" extra_about_right = %i\n", nextra); } } printf("ZZZ\n"); if (fPriority[2]==2) { for (id=0; id<fIbz; id++) { printf("%15.10f\n", fZb[id]); if (id == (fIbz-1)) break; printf("slice %i : %i\n", id, fIndZ[fOBz[id]]); if (fIndZ[fOBz[id]]) { slice = (UChar_t*)(&fIndZ[fOBz[id]+1]); for (i=0; i<nbytes; i++) { byte = slice[i]; for (bit=0; bit<8; bit++) { if (byte & (1<<bit)) printf(" %i ", 8*i+bit); } } printf("\n"); } GetExtraZ(id,kTRUE,nextra); printf(" extra_about_left = %i\n", nextra); GetExtraZ(id,kFALSE,nextra); printf(" extra_about_right = %i\n", nextra); } } } //----------------------------------------------------------------------------- void TGeoVoxelFinder::PrintVoxelLimits(Double_t *point) const { // print the voxel containing point Int_t im=0; if (fPriority[0]) { im = TMath::BinarySearch(fIbx, fXb, point[0]); if ((im==-1) || (im==fIbx-1)) { printf("Voxel X limits: OUT\n"); } else { printf("Voxel X limits: %g %g\n", fXb[im], fXb[im+1]); } } if (fPriority[1]) { im = TMath::BinarySearch(fIby, fYb, point[1]); if ((im==-1) || (im==fIby-1)) { printf("Voxel Y limits: OUT\n"); } else { printf("Voxel Y limits: %g %g\n", fYb[im], fYb[im+1]); } } if (fPriority[2]) { im = TMath::BinarySearch(fIbz, fZb, point[2]); if ((im==-1) || (im==fIbz-1)) { printf("Voxel Z limits: OUT\n"); } else { printf("Voxel Z limits: %g %g\n", fZb[im], fZb[im+1]); } } } //----------------------------------------------------------------------------- void TGeoVoxelFinder::Voxelize(Option_t * /*option*/) { // Voxelize attached volume according to option BuildVoxelLimits(); SortAll(); } //----------------------------------------------------------------------------- ClassImp(TGeoCylVoxels) //----------------------------------------------------------------------------- TGeoCylVoxels::TGeoCylVoxels() { // Default constructor } //----------------------------------------------------------------------------- TGeoCylVoxels::TGeoCylVoxels(TGeoVolume *vol) :TGeoVoxelFinder(vol) { // Constructor fVolume->SetCylVoxels(kTRUE); } //----------------------------------------------------------------------------- TGeoCylVoxels::~TGeoCylVoxels() { // Destructor } //----------------------------------------------------------------------------- void TGeoCylVoxels::BuildVoxelLimits() { //--- Compute boundary limits in R, Phi and Z coordinates for all daughters // of fVolume Int_t id; Int_t nd = fVolume->GetNdaughters(); TGeoNode *node; Double_t bcyl[4]; if (fBoxes) delete [] fBoxes; fNboxes = 6*nd; fBoxes = new Double_t[fNboxes]; if (fCheckList) delete [] fCheckList; fCheckList = new Int_t[nd]; Double_t vert[24]; Double_t pt[3]; Double_t xyz[6]; const Double_t *translation; TGeoBBox *box = 0; // Double_t *origin; Double_t dx, dy, dz, x0, y0; Double_t orig[3]; TGeoShape *shape = 0; TGeoMatrix *matrix = 0; // loop all daughters for (id=0; id<nd; id++) { node = fVolume->GetNode(id); // printf(" --%s--\n", node->GetName()); shape = node->GetVolume()->GetShape(); // shape->InspectShape(); box = (TGeoBBox*)shape; // origin = box->GetOrigin(); matrix = node->GetMatrix(); box->SetBoxPoints(&vert[0]); for (Int_t point=0; point<8; point++) { matrix->LocalToMaster(&vert[3*point], &pt[0]); if (!point) { xyz[0] = xyz[1] = pt[0]; xyz[2] = xyz[3] = pt[1]; xyz[4] = xyz[5] = pt[2]; continue; } for (Int_t j=0; j<3; j++) { if (pt[j] < xyz[2*j]) xyz[2*j]=pt[j]; if (pt[j] > xyz[2*j+1]) xyz[2*j+1]=pt[j]; } } dx = 0.5*(xyz[1]-xyz[0]); dy = 0.5*(xyz[3]-xyz[2]); dz = 0.5*(xyz[5]-xyz[4]); x0 = 0.5*(xyz[0]+xyz[1]); y0 = 0.5*(xyz[2]+xyz[3]); orig[0] = TMath::Abs(x0); orig[1] = TMath::Abs(y0); orig[2] = 0.5*(xyz[4]+xyz[5]); fBoxes[6*id+4] = orig[2]-dz; fBoxes[6*id+5] = orig[2]+dz; if (matrix->IsIdentity()) { // node has no rotation // printf(" identity\n"); shape->GetBoundingCylinder(&bcyl[0]); memcpy(fBoxes+6*id, &bcyl[0], 4*sizeof(Double_t)); } else { // matrix->Print(); if (matrix->IsRotAboutZ()) { // printf(" rotz\n"); // no rotation about other axis than Z translation = matrix->GetTranslation(); if ((TMath::Abs(translation[0])<1E-10) && (TMath::Abs(translation[1])<1E-10)) { // there is just a translation on Z (and possibly a rot. about Z) // printf(" Z transl.\n"); shape->GetBoundingCylinder(&bcyl[0]); // check if any rotation about Z if (matrix->IsRotation()) { // printf(" + rot\n"); // find phi rotation if ((bcyl[3]-bcyl[2])!=360) { Double_t phi = ((TGeoRotation *)matrix)->GetPhiRotation(); bcyl[2] += phi; bcyl[3] += phi; if (bcyl[2]<0) { bcyl[2] += 360.; bcyl[3] += 360.; } else { if (bcyl[2]>360.) { bcyl[2] -= 360.; bcyl[3] -=360.; } } } } } else { // translation is other than Z // printf(" gen. translation\n"); memset(&xyz[0], 0, 6*sizeof(Double_t)); // origin of mother to local frame matrix->MasterToLocal(&xyz[0], &xyz[3]); x0 = TMath::Abs(xyz[3]); y0 = TMath::Abs(xyz[4]); dx = box->GetDX(); dy = box->GetDY(); bcyl[1] = (x0+dx)*(x0+dx)+(y0+dy)*(y0+dy); if (x0<dx) { // origin in X range // printf(" inside X\n"); if (y0<dy) { // origin also in Y range bcyl[0] = 0.; bcyl[2] = 0.; bcyl[3] = 360.; } else { // origin outside Y range // printf(" outside Y\n"); bcyl[0] = y0-dy; bcyl[0] *= bcyl[0]; // convert phi limits to MARS if (xyz[4]>0) { xyz[3] = -dx; xyz[4] = dy; } else { xyz[3] = dx; xyz[4] = -dy; } matrix->LocalToMaster(&xyz[3], &xyz[0]); // printf(" at phi1: %g %g\n", xyz[0], xyz[1]); bcyl[2] = TMath::ATan2(xyz[1], xyz[0])*TMath::RadToDeg(); xyz[3] = -xyz[3]; matrix->LocalToMaster(&xyz[3], &xyz[0]); // printf(" at phi2: %g %g\n", xyz[0], xyz[1]); bcyl[3] = TMath::ATan2(xyz[1], xyz[0])*TMath::RadToDeg(); if (bcyl[2]<0) bcyl[2]+=360.; while (bcyl[3]<bcyl[2]) bcyl[3]+=360.; } } else { // origin outside X range // printf(" outside X\n"); if (y0<dy) { // origin in Y range // printf(" inside Y\n"); bcyl[0] = x0-dx; bcyl[0] *= bcyl[0]; // convert phi limits to MARS if (xyz[3]>0) { xyz[3] = dx; xyz[4] = dy; } else { xyz[3] = -dx; xyz[4] = -dy; } matrix->LocalToMaster(&xyz[3], &xyz[0]); // printf(" at phi1: %g %g\n", xyz[0], xyz[1]); bcyl[2] = TMath::ATan2(xyz[1], xyz[0])*TMath::RadToDeg(); xyz[4] = -xyz[4]; matrix->LocalToMaster(&xyz[3], &xyz[0]); // printf(" at phi2: %g %g\n", xyz[0], xyz[1]); bcyl[3] = TMath::ATan2(xyz[1], xyz[0])*TMath::RadToDeg(); if (bcyl[2]<0) bcyl[2]+=360.; while (bcyl[3]<bcyl[2]) bcyl[3]+=360.; } else { // origin outside both X and Y range // printf(" outside XY\n"); bcyl[0] = (x0-dx)*(x0-dx)+(y0-dy)*(y0-dy); // convert phi limits to MARS if (xyz[3]>0) { xyz[3] = xyz[4]; xyz[4] = dy; xyz[3]=(xyz[3]>0)?-dx:dx; } else { xyz[3] = xyz[4]; xyz[4] = -dy; xyz[3]=(xyz[3]>0)?-dx:dx; } matrix->LocalToMaster(&xyz[3], &xyz[0]); // printf(" at phi1: %g %g\n", xyz[0], xyz[1]); bcyl[2] = TMath::ATan2(xyz[1], xyz[0])*TMath::RadToDeg(); xyz[3] = -xyz[3]; xyz[4] = -xyz[4]; matrix->LocalToMaster(&xyz[3], &xyz[0]); // printf(" at phi2: %g %g\n", xyz[0], xyz[1]); bcyl[3] = TMath::ATan2(xyz[1], xyz[0])*TMath::RadToDeg(); if (bcyl[2]<0) bcyl[2]+=360.; while (bcyl[3]<bcyl[2]) bcyl[3]+=360.; } } } // printf(" ---copy param\n"); memcpy(&fBoxes[6*id], &bcyl[0], 4*sizeof(Double_t)); } else { // general rotation matrix (not only about Z) // printf("General tranformation for node %s\n", node->GetName()); bcyl[1] = (orig[0]+dx)*(orig[0]+dx)+(orig[1]+dy)*(orig[1]+dy); if (orig[0]<dx) { if (orig[1]<dy) { bcyl[0] = 0.; bcyl[2] = 0.; bcyl[3] = 360.; } else { bcyl[0] = orig[1]-dy; bcyl[0] *= bcyl[0]; if (y0>0) { bcyl[2] = TMath::RadToDeg()*TMath::ATan2(xyz[2], xyz[1]); bcyl[3] = TMath::RadToDeg()*TMath::ATan2(xyz[2], xyz[0]); } else { bcyl[2] = TMath::RadToDeg()*TMath::ATan2(xyz[3], xyz[0]); bcyl[3] = TMath::RadToDeg()*TMath::ATan2(xyz[3], xyz[1]); } if (bcyl[2]<0) bcyl[2]+=360.; while (bcyl[3]<bcyl[2]) bcyl[3]+=360.; } } else { if (orig[1]<dy) { bcyl[0] = orig[0]-dx; bcyl[0] *= bcyl[0]; if (x0>0) { bcyl[2] = TMath::RadToDeg()*TMath::ATan2(xyz[2], xyz[0]); bcyl[3] = TMath::RadToDeg()*TMath::ATan2(xyz[3], xyz[0]); } else { bcyl[2] = TMath::RadToDeg()*TMath::ATan2(xyz[3], xyz[1]); bcyl[3] = TMath::RadToDeg()*TMath::ATan2(xyz[2], xyz[1]); } if (bcyl[2]<0) bcyl[2]+=360.; while (bcyl[3]<bcyl[2]) bcyl[3]+=360.; } else { bcyl[0] = (orig[0]-dx)*(orig[0]-dx)+(orig[1]-dy)*(orig[1]-dy); Int_t indx, indy; indy = (x0>0)?0:1; indx = (y0>0)?1:0; bcyl[2] = TMath::RadToDeg()*TMath::ATan2(xyz[indy+2], xyz[indx]); bcyl[3] = TMath::RadToDeg()*TMath::ATan2(xyz[3-indy], xyz[1-indx]); if (bcyl[2]<0) bcyl[2]+=360.; while (bcyl[3]<bcyl[2]) bcyl[3]+=360.; } } } memcpy(&fBoxes[6*id], &bcyl[0], 4*sizeof(Double_t)); } // printf("Limits for %s\n", node->GetName()); // printf(" R : %g %g\n", fBoxes[6*id], fBoxes[6*id+1]); // printf(" Phi : %g %g\n", fBoxes[6*id+2], fBoxes[6*id+3]); // printf(" Z : %g %g\n", fBoxes[6*id+4], fBoxes[6*id+5]); } } //----------------------------------------------------------------------------- Double_t TGeoCylVoxels::Efficiency() { //--- Compute voxelization efficiency. return 0; } //----------------------------------------------------------------------------- void TGeoCylVoxels::FindOverlaps(Int_t inode) const { // create the list of nodes for which the bboxes overlap with inode's bbox if (!fBoxes) return; Double_t xmin, xmax, ymin, ymax, zmin, zmax; Double_t xmin1, xmax1, ymin1, ymax1, zmin1, zmax1; Double_t ddx1, ddx2; Int_t nd = fVolume->GetNdaughters(); Int_t *ovlps = 0; Int_t *otmp = new Int_t[nd-1]; Int_t novlp = 0; TGeoNode *node = fVolume->GetNode(inode); xmin = fBoxes[6*inode]; xmax = fBoxes[6*inode+1]; ymin = fBoxes[6*inode+2]; ymax = fBoxes[6*inode+3]; zmin = fBoxes[6*inode+4]; zmax = fBoxes[6*inode+5]; // printf("overlaps for MANY node %s\n", node->GetName()); // printf("xmin=%g xmax=%g\n", xmin, xmax); // printf("ymin=%g ymax=%g\n", ymin, ymax); // printf("zmin=%g zmax=%g\n", zmin, zmax); Bool_t in = kFALSE; //TGeoNode *node1; // loop on brothers for (Int_t ib=0; ib<nd; ib++) { if (ib == inode) continue; // everyone overlaps with itself in = kFALSE; //node1 = fVolume->GetNode(ib); xmin1 = fBoxes[6*ib]; xmax1 = fBoxes[6*ib+1]; ymin1 = fBoxes[6*ib+2]; ymax1 = fBoxes[6*ib+3]; zmin1 = fBoxes[6*ib+4]; zmax1 = fBoxes[6*ib+5]; // printf(" node %s\n", node1->GetName()); // printf(" xmin1=%g xmax1=%g\n", xmin1, xmax1); // printf(" ymin1=%g ymax1=%g\n", ymin1, ymax1); // printf(" zmin1=%g zmax1=%g\n", zmin1, zmax1); ddx1 = TMath::Abs(xmin1-xmax); ddx2 = TMath::Abs(xmax1-xmin); if ((ddx1<1E-6)||(ddx2<1E-6)) continue; // if ((xmin1==xmin)||(xmax1==xmax)) in = kTRUE; if (((xmin1>xmin)&&(xmin1<xmax))||((xmax1>xmin)&&(xmax1<xmax)) || ((xmin>xmin1)&&(xmin<xmax1))||((xmax>xmin1)&&(xmax<xmax1))) in = kTRUE; if (!in) continue; // printf("x overlap...\n"); in = kFALSE; if (ymax<360.) { in = (IntersectIntervals(ymin, ymax, ymin1, ymax1)>0)?kTRUE:kFALSE; } else { if (ymax1<360.) { in = (IntersectIntervals(ymin1, ymax1, ymin, ymax)>0)?kTRUE:kFALSE; } else { in = (IntersectIntervals(ymin1, 360., ymin, ymax)>0)?kTRUE:kFALSE; if (!in) in = (IntersectIntervals(360., ymax1, ymin, ymax)>0)?kTRUE:kFALSE; } } if (!in) continue; // printf("y overlap...\n"); in = kFALSE; ddx1 = TMath::Abs(zmin1-zmax); ddx2 = TMath::Abs(zmax1-zmin); if ((ddx1<1E-12)||(ddx2<1E-12)) continue; // if ((zmin1==zmin)||(zmax1==zmax)) in = kTRUE; if (((zmin1>zmin)&&(zmin1<zmax))||((zmax1>zmin)&&(zmax1<zmax)) || ((zmin>zmin1)&&(zmin<zmax1))||((zmax>zmin1)&&(zmax<zmax1))) in = kTRUE; if (!in) continue; // printf("Overlapping %i\n", ib); otmp[novlp++] = ib; } if (!novlp) { // printf("---no overlaps for MANY node %s\n", node->GetName()); node->SetOverlaps(ovlps, 1); return; } ovlps = new Int_t[novlp]; memcpy(ovlps, otmp, novlp*sizeof(Int_t)); delete [] otmp; node->SetOverlaps(ovlps, novlp); // printf("Overlaps for MANY node %s : %i\n", node->GetName(), novlp); } //----------------------------------------------------------------------------- Int_t *TGeoCylVoxels::GetCheckList(Double_t *point, Int_t &nelem) { //--- Get the list of nodes possibly containing a given point. // convert the point to cylindrical coordinates Double_t ptcyl[3]; ptcyl[0] = point[0]*point[0]+point[1]*point[1]; if (fPriority[1]) { ptcyl[1] = TMath::ATan2(point[1], point[0])*TMath::RadToDeg(); if (ptcyl[1]<0) ptcyl[1]+=360.; } ptcyl[2] = point[2]; return TGeoVoxelFinder::GetCheckList(&ptcyl[0], nelem); } //----------------------------------------------------------------------------- Int_t *TGeoCylVoxels::GetNextVoxel(Double_t * /*point*/, Double_t * /*dir*/, Int_t & /*ncheck*/) { //--- Get the list of nodes possibly crossed by a given ray. return 0; } //----------------------------------------------------------------------------- Int_t TGeoCylVoxels::IntersectIntervals(Double_t vox1, Double_t vox2, Double_t phi1, Double_t phi2) const { // Intersect a phi voxel interval in range (0,360) with a node phi interval of // extended range. Returns 0 if no intersection, 1 if they do intersect and // 2 if they are identical. if ((vox2-vox1)==360.) { if ((phi2-phi1)==360.) return 2; return 1; } Double_t d11, d12; // check if first phi limit correspond to voxel d11 = phi1-vox1; if (TMath::Abs(d11)<1E-8) d11=0; d12 = phi1-vox2; if (TMath::Abs(d12)<1E-8) d12=0; // check if second phi limit is in range (0, 360) if (phi2>360.) { if (d11>=0) { if (d12<0) return 1; d11 = vox1-phi2+360.; if (d11>=0) return 0; return 1; } return 1; } // both intervals are in range (0, 360) Double_t d22 = phi2-vox2; if (TMath::Abs(d22)<1E-8) d22=0; if ((d11==0.) && (d22==0.)) return 2; if (d11>=0) { if (d12<0) return 1; return 0; } Double_t d21 = phi2-vox1; if (TMath::Abs(d21)<1E-8) d21=0; if (d21>0) return 1; return 0; } //----------------------------------------------------------------------------- void TGeoCylVoxels::Print(Option_t *) const { // Print info about voxels. Int_t id; printf("Voxels for volume %s (nd=%i)\n", fVolume->GetName(), fVolume->GetNdaughters()); printf("priority : x=%i y=%i z=%i\n", fPriority[0], fPriority[1], fPriority[2]); // return; printf("--- R voxels ---\n"); if (fPriority[0]) { for (id=0; id<fIbx; id++) { printf("%15.10f\n",TMath::Sqrt(fXb[id])); if (id == (fIbx-1)) break; printf("slice %i : %i\n", id, fIndX[fOBx[id]]); for (Int_t j=0;j<fIndX[fOBx[id]]; j++) { printf("%s low, high: %15.10f --- %15.10f \n", fVolume->GetNode(fIndX[fOBx[id]+j+1])->GetName(), TMath::Sqrt(fBoxes[6*fIndX[fOBx[id]+j+1]]), TMath::Sqrt(fBoxes[6*fIndX[fOBx[id]+j+1]+1])); } } } printf("--- Phi voxels ---\n"); if (fPriority[1]) { for (id=0; id<fIby; id++) { printf("%15.10f\n", fYb[id]); if (id == (fIby-1)) break; printf("slice %i : %i\n", id, fIndY[fOBy[id]]); for (Int_t j=0;j<fIndY[fOBy[id]]; j++) { printf("%s low, high: %15.10f --- %15.10f \n", fVolume->GetNode(fIndY[fOBy[id]+j+1])->GetName(), fBoxes[6*fIndY[fOBy[id]+j+1]+2], fBoxes[6*fIndY[fOBy[id]+j+1]+3]); } } } printf(" ---Z voxels---\n"); if (fPriority[2]) { for (id=0; id<fIbz; id++) { printf("%15.10f\n", fZb[id]); if (id == (fIbz-1)) break; printf("slice %i : %i\n", id, fIndZ[fOBz[id]]); for (Int_t j=0;j<fIndZ[fOBz[id]]; j++) { printf("%s low, high: %15.10f --- %15.10f \n", fVolume->GetNode(fIndZ[fOBz[id]+j+1])->GetName(), fBoxes[6*fIndZ[fOBz[id]+j+1]+4], fBoxes[6*fIndZ[fOBz[id]+j+1]+5]); } } } } //----------------------------------------------------------------------------- void TGeoCylVoxels::SortAll(Option_t *) { // Order voxels along R, Phi, Z. // printf("Sorting voxels for %s\n", fVolume->GetName()); Int_t nd = fVolume->GetNdaughters(); if (!nd) return; Double_t *boundaries = new Double_t[6*nd]; Double_t rmin, rmax, pmin, pmax, zmin, zmax; Double_t phi2; Double_t bcyl[4]; TGeoShape *shape = fVolume->GetShape(); shape->GetBoundingCylinder(&bcyl[0]); TGeoBBox *box = (TGeoBBox*)fVolume->GetShape(); // compute ranges on R, Phi, Z according to volume cylindrical limits rmin = bcyl[0]; rmax = bcyl[1]; pmin = bcyl[2]; pmax = bcyl[3]; zmin = (box->GetOrigin())[2] - box->GetDZ(); zmax = (box->GetOrigin())[2] + box->GetDZ(); if ((rmin>=rmax) || (pmin>=pmax) || (zmin>=zmax)) { Error("SortAll", "wrong bounding cylinder"); printf("### volume was : %s\n", fVolume->GetName()); SetInvalid(); return; } Int_t id; // compute boundaries coordinates on R, Phi and Z for (id=0; id<nd; id++) { // r boundaries boundaries[2*id] = fBoxes[6*id]; boundaries[2*id+1] = fBoxes[6*id+1]; // phi boundaries boundaries[2*id+2*nd] = fBoxes[6*id+2]; phi2 = fBoxes[6*id+3]; if (phi2>360.) phi2-=360.; boundaries[2*id+2*nd+1] = phi2; // z boundaries boundaries[2*id+4*nd] = fBoxes[6*id+4]; boundaries[2*id+4*nd+1] = fBoxes[6*id+5]; } Int_t *index = new Int_t[2*nd]; Int_t *ind = new Int_t[(2*nd+2)*(nd+1)]; // ind[fOBx[i]] = ndghts in slice fInd[i]--fInd[i+1] Double_t *temp = new Double_t[2*nd+2]; Int_t current = 0; Double_t xxmin, xxmax, xbmin, xbmax, ddx1, ddx2; Int_t last; Double_t db, dbv; // sort r boundaries Int_t ib = 0; TMath::Sort(2*nd, &boundaries[0], &index[0], kFALSE); last = index[0]; db = boundaries[last+1-2*last%2]-boundaries[last]; temp[ib++] = boundaries[last]; // compact common boundaries for (id=1; id<2*nd; id++) { dbv = boundaries[index[id]]-boundaries[last]; last = index[id]; if (dbv>1E-6) { // we have to generate a new boundary temp[ib++] = boundaries[last]; db = boundaries[last+1-2*last%2]-boundaries[last]; } else { if (db<0) { // just ignore this boundary that have to be compacted with an other temp[ib-1] = boundaries[last]; } } } // now find priority if (ib < 2) { Error("SortAll", "less than 2 boundaries on R !"); printf("### volume was : %s\n", fVolume->GetName()); SetInvalid(); return; } if (ib == 2) { // check range if (((temp[0]-rmin)<1E-8) && ((temp[1]-rmax)>-1E-8)) { // ordering on this axis makes no sense. Clear all arrays. fPriority[0] = 0; if (fIndX) delete [] fIndX; fIndX = 0; fNx = 0; if (fXb) delete [] fXb; fXb = 0; fIbx = 0; if (fOBx) delete [] fOBx; fOBx = 0; fNox = 0; } else { fPriority[0] = 1; // all in one slice } } else { fPriority[0] = 2; // check all } // store compacted boundaries if (fPriority[0]) { if (fXb) delete [] fXb; fXb = new Double_t[ib]; memcpy(fXb, &temp[0], ib*sizeof(Double_t)); fIbx = ib; //now build the lists of nodes in each slice memset(ind, 0, (2*nd+2)*(nd+1)*sizeof(Int_t)); // ind[fOBx[i]+k] = index of dght k (k<ndghts) if (fOBx) delete [] fOBx; fNox = fIbx-1; fOBx = new Int_t[fNox]; // offsets in ind for (id=0; id<fNox; id++) { fOBx[id] = current; // offset of dght list ind[current] = 0; // ndght in this slice xxmin = fXb[id]; xxmax = fXb[id+1]; for (Int_t ic=0; ic<nd; ic++) { xbmin = fBoxes[6*ic]; xbmax = fBoxes[6*ic+1]; ddx1 = TMath::Abs(xbmin-xxmax); ddx2 = TMath::Abs(xbmax-xxmin); if ((xbmin==xxmin)||(xbmax==xxmax)) { ind[current]++; ind[current+ind[current]] = ic; continue; } if ((ddx1<1E-8)||(ddx2<1E-8)) continue; if (((xbmin>xxmin)&&(xbmin<xxmax))||((xbmax>xxmin)&&(xbmax<xxmax)) || ((xxmin>xbmin)&&(xxmin<xbmax))||((xxmax>xbmin)&&(xxmax<xbmax))) { // daughter ic in interval ind[current]++; ind[current+ind[current]] = ic; } } current += ind[current]+1; } if (fIndX) delete [] fIndX; fNx = current; fIndX = new Int_t[current]; memcpy(fIndX, &ind[0], current*sizeof(Int_t)); } // sort Phi boundaries Int_t intersect; ib = 0; TMath::Sort(2*nd, &boundaries[2*nd], &index[0], kFALSE); temp[ib++] = 0.; // always store phi=0 and phi=360 boundaries // compact common boundaries for (id=0; id<2*nd; id++) { if (TMath::Abs(temp[ib-1]-boundaries[2*nd+index[id]])>1E-8) temp[ib++]=boundaries[2*nd+index[id]]; } if (temp[ib-1]!=360.) temp[ib++]=360.; // now find priority on Phi if (ib < 2) { Error("SortAll", "less than 2 boundaries on Phi !"); printf("### volume was : %s\n", fVolume->GetName()); SetInvalid(); return; } if (ib == 2) { // check range intersect = IntersectIntervals(temp[0], temp[1], pmin, pmax); if (intersect==2) { // ordering on this axis makes no sense. Clear all arrays. fPriority[1] = 0; if (fIndY) delete [] fIndY; fIndY = 0; fNy = 0; if (fYb) delete [] fYb; fYb = 0; fIby = 0; if (fOBy) delete [] fOBy; fOBy = 0; fNoy = 0; } else { fPriority[1] = 1; // all in one slice } } else { fPriority[1] = 2; // check all } if (fPriority[1]) { // store compacted boundaries if (fYb) delete [] fYb; fYb = new Double_t[ib]; memcpy(fYb, &temp[0], ib*sizeof(Double_t)); fIby = ib; memset(ind, 0, (2*nd+2)*(nd+1)*sizeof(Int_t)); current = 0; if (fOBy) delete [] fOBy; fNoy = fIby-1; fOBy = new Int_t[fNoy]; // offsets in ind for (id=0; id<fNoy; id++) { // loop voxels fOBy[id] = current; // offset of dght list ind[current] = 0; // ndght in this slice xxmin = fYb[id]; // voxel limits xxmax = fYb[id+1]; for (Int_t ic=0; ic<nd; ic++) { // loop daughters xbmin = fBoxes[6*ic+2]; xbmax = fBoxes[6*ic+3]; intersect = IntersectIntervals(xxmin, xxmax, xbmin, xbmax); if (intersect) { ind[current]++; ind[current+ind[current]] = ic; } } current += ind[current]+1; } if (fIndY) delete [] fIndY; fNy = current; fIndY = new Int_t[current]; memcpy(fIndY, &ind[0], current*sizeof(Int_t)); } // sort z boundaries ib = 0; TMath::Sort(2*nd, &boundaries[4*nd], &index[0], kFALSE); temp[ib++] = boundaries[4*nd+index[0]]; // compact common boundaries for (id=1; id<2*nd; id++) { if ((TMath::Abs(temp[ib-1]-boundaries[4*nd+index[id]]))>1E-10) temp[ib++]=boundaries[4*nd+index[id]]; } // now find priority on Z if (ib < 2) { Error("SortAll", "less than 2 boundaries on Z !"); printf("### volume was : %s\n", fVolume->GetName()); SetInvalid(); return; } if (ib == 2) { // check range if (((temp[0]-zmin)<1E-10) && ((temp[1]-zmax)>-1E-10)) { // ordering on this axis makes no sense. Clear all arrays. fPriority[2] = 0; if (fIndZ) delete [] fIndZ; fIndZ = 0; fNz = 0; if (fZb) delete [] fZb; fZb = 0; fIbz = 0; if (fOBz) delete [] fOBz; fOBz = 0; fNoz = 0; } else { fPriority[2] = 1; // all in one slice } } else { fPriority[2] = 2; // check all } if (fPriority[2]) { // store compacted boundaries if (fZb) delete [] fZb; fZb = new Double_t[ib]; memcpy(fZb, &temp[0], (ib)*sizeof(Double_t)); fIbz = ib; memset(ind, 0, (2*nd+2)*(nd+1)*sizeof(Int_t)); current = 0; if (fOBz) delete [] fOBz; fNoz = fIbz-1; fOBz = new Int_t[fNoz]; // offsets in ind for (id=0; id<fNoz; id++) { fOBz[id] = current; // offset of dght list ind[current] = 0; // ndght in this slice xxmin = fZb[id]; xxmax = fZb[id+1]; for (Int_t ic=0; ic<nd; ic++) { xbmin = fBoxes[6*ic+4]; xbmax = fBoxes[6*ic+5]; ddx1 = TMath::Abs(xbmin-xxmax); ddx2 = TMath::Abs(xbmax-xxmin); if ((xbmin==xxmin)||(xbmax==xxmax)) { ind[current]++; ind[current+ind[current]] = ic; continue; } if ((ddx1<1E-12)||(ddx2<1E-12)) continue; if (((xbmin>xxmin)&&(xbmin<xxmax))||((xbmax>xxmin)&&(xbmax<xxmax)) || ((xxmin>xbmin)&&(xxmin<xbmax))||((xxmax>xbmin)&&(xxmax<xbmax))) { // daughter ic in interval ind[current]++; ind[current+ind[current]] = ic; } } current += ind[current]+1; } if (fIndZ) delete [] fIndZ; fNz = current; fIndZ = new Int_t[current]; memcpy(fIndZ, &ind[0], current*sizeof(Int_t)); } delete [] boundaries; boundaries=0; delete [] index; index=0; delete [] temp; temp=0; delete [] ind; // Print(); if ((!fPriority[0]) && (!fPriority[1]) && (!fPriority[2])) { fVolume->SetVoxelFinder(0); delete this; } else { // Efficiency(); } } //----------------------------------------------------------------------------- void TGeoCylVoxels::Voxelize(Option_t *) { //--- Voxelize fVolume. // printf("Voxelizing %s\n", fVolume->GetName()); BuildVoxelLimits(); SortAll(); } ClassImp(TGeoFullVoxels) //----------------------------------------------------------------------------- TGeoFullVoxels::TGeoFullVoxels() { // Default constructor fNvoxels = 0; fNvx = 0; fNvy = 0; fNvz = 0; fVox = 0; } //----------------------------------------------------------------------------- TGeoFullVoxels::TGeoFullVoxels(TGeoVolume *vol) :TGeoVoxelFinder(vol) { // Constructor fNvoxels = 0; fNvx = 0; fNvy = 0; fNvz = 0; fVox = 0; } //----------------------------------------------------------------------------- TGeoFullVoxels::~TGeoFullVoxels() { // Destructor if (fVox) delete [] fVox; } //----------------------------------------------------------------------------- void TGeoFullVoxels::Voxelize(Option_t *) { //--- Voxelize fVolume. // printf("Voxelizing %s\n", fVolume->GetName()); TGeoVoxelFinder::Voxelize(); if (IsInvalid()) return; fNvx = fNvy = fNvz = 1; if (fPriority[0]) fNvx = fIbx-1; if (fPriority[1]) fNvy = fIby-1; if (fPriority[2]) fNvz = fIbz-1; fNvoxels = fNvx*fNvy*fNvz; if (fNvoxels <= 0) { SetInvalid(); return; } fVox = new UChar_t[fNvoxels]; // Intersect slices to get candidates Int_t ptrx=0, ptry=0, ptrz=0; // index of voxel (i,j,k) UChar_t *slice1 = 0; UChar_t *slice2 = 0; UChar_t *slice3 = 0; Int_t i,j,k; if (fPriority[0]==2) { for (i=0; i<fNvx; i++) { ptrx = i*fNvy*fNvz; slice1 = (UChar_t*)(&fIndX[fOBx[i]+1]); if (fPriority[1]==2) { for (j=0; j<fNvy; j++) { ptry = ptrx + j*fNvz; slice2 = (UChar_t*)(&fIndY[fOBy[j]+1]); if (fPriority[2]==2) { for (k=0; k<fNvz; k++) { ptrz = ptry + k; slice3 = (UChar_t*)(&fIndZ[fOBz[k]+1]); fVox[ptrz] = slice1[0] & slice2[0] & slice3[0]; } } else { fVox[ptry] = slice1[0] & slice2[0]; } } } else { if (fPriority[2]==2) { for (k=0; k<fNvz; k++) { ptrz = ptrx + k; slice3 = (UChar_t*)(&fIndZ[fOBz[k]+1]); fVox[ptrz] = slice1[0] & slice3[0]; } } else { fVox[ptrx] = slice1[0]; } } } } else { if (fPriority[1]==2) { for (j=0; j<fNvy; j++) { ptry = j*fNvz; slice2 = (UChar_t*)(&fIndY[fOBy[j]+1]); if (fPriority[2]==2) { for (k=0; k<fNvz; k++) { ptrz = ptry + k; slice3 = (UChar_t*)(&fIndZ[fOBz[k]+1]); fVox[ptrz] = slice2[0] & slice3[0]; } } else { fVox[ptry] = slice2[0]; } } } else { if (fPriority[2]==2) { for (k=0; k<fNvz; k++) { ptrz = k; slice3 = (UChar_t*)(&fIndZ[fOBz[k]+1]); fVox[ptrz] = slice3[0]; } } } } } //----------------------------------------------------------------------------- Int_t *TGeoFullVoxels::GetVoxelCandidates(Int_t i, Int_t j, Int_t k, Int_t &ncheck) { // get the list of candidates in voxel (i,j,k) - no check ncheck = 0; Int_t nd = fVolume->GetNdaughters(); UChar_t *vox = GetVoxel(i,j,k); UChar_t byte = vox[0]; if (!vox[0]) return 0; for (Int_t bit=0; bit<nd; bit++) { if (byte & (1<<bit)) fCheckList[ncheck++] = bit; } return fCheckList; } //----------------------------------------------------------------------------- Int_t *TGeoFullVoxels::GetCheckList(Double_t *point, Int_t &nelem) { // get the list of daughter indices for which point is inside their bbox nelem = fNcandidates = 0; Int_t im; UChar_t *slice; UChar_t byte = 0xFF; if (fPriority[2]) { im = TMath::BinarySearch(fIbz, fZb, point[2]); if ((im==-1) || (im==fIbz-1)) return 0; if (fPriority[2]==2) { slice = (UChar_t*)(&fIndZ[fOBz[im]+1]); byte &= slice[0]; if (!byte) return 0; } } if (fPriority[0]) { im = TMath::BinarySearch(fIbx, fXb, point[0]); if ((im==-1) || (im==fIbx-1)) return 0; if (fPriority[0]==2) { slice = (UChar_t*)(&fIndX[fOBx[im]+1]); byte &= slice[0]; if (!byte) return 0; } } if (fPriority[1]) { im = TMath::BinarySearch(fIby, fYb, point[1]); if ((im==-1) || (im==fIby-1)) return 0; if (fPriority[1]==2) { slice = (UChar_t*)(&fIndY[fOBy[im]+1]); byte &= slice[0]; if (!byte) return 0; } } Int_t nd = fVolume->GetNdaughters(); for (Int_t i=0; i<nd; i++) { if (byte & (1<<i)) fCheckList[fNcandidates++] = i; } nelem = fNcandidates; return fCheckList; } //----------------------------------------------------------------------------- void TGeoFullVoxels::Print(Option_t *) const { }