From 7672e15e77ce6b71fa0c532a56e1699050833e04 Mon Sep 17 00:00:00 2001
From: Rene Brun <Rene.Brun@cern.ch>
Date: Tue, 20 Jan 2004 15:44:33 +0000
Subject: [PATCH] Several improvements from Andrei GHeata.
additional option "i" (reset by default) to CloseGeometry() to have the
node id array built or not.
git-svn-id: http://root.cern.ch/svn/root/trunk@7952 27541ba8-7e3a-0410-8455-c3a389f83636
---
geom/Module.mk | 3 +-
geom/inc/LinkDef1.h | 4 +-
geom/inc/TGeoCache.h | 19 ++-
geom/inc/TGeoManager.h | 9 +-
geom/inc/TGeoMatrix.h | 11 +-
geom/inc/TGeoVoxelFinder.h | 43 +++++-
geom/src/TGeoCache.cxx | 126 +++++++++++-----
geom/src/TGeoManager.cxx | 48 ++++--
geom/src/TGeoTube.cxx | 4 +-
geom/src/TGeoVolume.cxx | 21 ++-
geom/src/TGeoVoxelFinder.cxx | 282 +++++++++++++++++++++++++++++++++--
11 files changed, 475 insertions(+), 95 deletions(-)
diff --git a/geom/Module.mk b/geom/Module.mk
index 13d7be574d6..e7c47348ae5 100644
--- a/geom/Module.mk
+++ b/geom/Module.mk
@@ -30,7 +30,8 @@ GEOMH1 := TGeoAtt.h TGeoBoolNode.h \
TGeoEltu.h TGeoCone.h TGeoPcon.h \
TGeoPgon.h TGeoArb8.h TGeoTrd1.h TGeoTrd2.h \
TGeoManager.h TGeoCompositeShape.h \
- TVirtualGeoPainter.h TVirtualGeoTrack.h
+ TVirtualGeoPainter.h TVirtualGeoTrack.h \
+ TGeoPolygon.h
GEOMH2 := TGeoPatternFinder.h TGeoCache.h
GEOMH1 := $(patsubst %,$(MODDIRI)/%,$(GEOMH1))
GEOMH2 := $(patsubst %,$(MODDIRI)/%,$(GEOMH2))
diff --git a/geom/inc/LinkDef1.h b/geom/inc/LinkDef1.h
index cace0039ca8..e866a1cfe85 100644
--- a/geom/inc/LinkDef1.h
+++ b/geom/inc/LinkDef1.h
@@ -1,4 +1,4 @@
-// @(#)root/geom:$Name: $:$Id: LinkDef1.h,v 1.1 2003/09/09 16:56:48 brun Exp $
+// @(#)root/geom:$Name: $:$Id: LinkDef1.h,v 1.2 2003/10/20 08:46:33 brun Exp $
// Author : Andrei Gheata 10/06/02
/*************************************************************************
* Copyright (C) 1995-2000, Rene Brun and Fons Rademakers. *
@@ -35,6 +35,7 @@
#pragma link C++ class TGeoIdentity+;
#pragma link C++ class TGeoVoxelFinder+;
#pragma link C++ class TGeoCylVoxels+;
+#pragma link C++ class TGeoFullVoxels+;
#pragma link C++ class TGeoShape+;
#pragma link C++ class TGeoBBox+;
#pragma link C++ class TGeoPara+;
@@ -54,6 +55,7 @@
#pragma link C++ class TGeoTrd1+;
#pragma link C++ class TGeoTrd2+;
#pragma link C++ class TGeoCompositeShape+;
+#pragma link C++ class TGeoPolygon+;
#pragma link C++ class TGeoVolume-;
#pragma link C++ class TGeoVolumeAssembly+;
#pragma link C++ class TGeoVolumeMulti+;
diff --git a/geom/inc/TGeoCache.h b/geom/inc/TGeoCache.h
index ab3a0294394..1eb5d6cf22a 100644
--- a/geom/inc/TGeoCache.h
+++ b/geom/inc/TGeoCache.h
@@ -1,4 +1,4 @@
-// @(#)root/geom:$Name: $:$Id: TGeoCache.h,v 1.17 2003/10/01 17:53:11 brun Exp $
+// @(#)root/geom:$Name: $:$Id: TGeoCache.h,v 1.18 2003/11/11 15:44:28 brun Exp $
// Author: Andrei Gheata 18/03/02
/*************************************************************************
@@ -43,6 +43,7 @@ class TGeoCacheState : public TObject
{
protected:
Int_t fLevel; // level in the current branch
+ Int_t fStart; // start level
Int_t fIdBranch[30]; // ID branch
Double_t *fPoint; // last point in master frame
Bool_t fOverlapping; // overlap flag
@@ -55,7 +56,7 @@ public:
TGeoCacheState(Int_t capacity);
virtual ~TGeoCacheState();
- virtual void SetState(Int_t level, Bool_t ovlp, Double_t *point=0);
+ virtual void SetState(Int_t level, Int_t startlevel, Bool_t ovlp, Double_t *point=0);
virtual Bool_t GetState(Int_t &level, Double_t *point) const;
ClassDef(TGeoCacheState, 1) // class storing the cache state
@@ -78,7 +79,7 @@ public:
TGeoCacheStateDummy(Int_t capacity);
virtual ~TGeoCacheStateDummy();
- virtual void SetState(Int_t level, Bool_t ovlp, Double_t *point=0);
+ virtual void SetState(Int_t level, Int_t startlevel, Bool_t ovlp, Double_t *point=0);
virtual Bool_t GetState(Int_t &level, Double_t *point) const;
ClassDef(TGeoCacheStateDummy, 1) // class storing the cache state
@@ -129,7 +130,8 @@ protected:
public:
TGeoNodeCache();
- TGeoNodeCache(Int_t size);
+ TGeoNodeCache(Bool_t nodeid);
+ TGeoNodeCache(Int_t size, Bool_t nodeid=kFALSE);
virtual ~TGeoNodeCache();
Int_t AddNode(TGeoNode *node);
@@ -146,7 +148,7 @@ public:
virtual void ClearDaughter(Int_t index);
virtual void ClearNode(Int_t nindex);
virtual void Compact();
- void FillIdBranch(const Int_t *br) {memcpy(fIdBranch,br,(fLevel+1)*sizeof(Int_t)); fIndex=fIdBranch[fLevel];}
+ void FillIdBranch(const Int_t *br, Int_t startlevel=0) {memcpy(fIdBranch+startlevel,br,(fLevel+1-startlevel)*sizeof(Int_t)); fIndex=fIdBranch[fLevel];}
const Int_t *GetIdBranch() const {return fIdBranch;}
virtual void DeleteCaches();
virtual Bool_t DumpNodes();
@@ -157,7 +159,7 @@ public:
virtual void *GetMatrices() const {return fMatrices;}
virtual TGeoHMatrix *GetCurrentMatrix() const;
Int_t GetCurrentNode() const {return fCurrentNode;}
- Int_t GetCurrentNodeId() const {return (fNodeIdArray)?fNodeIdArray[fIndex]:-1;}
+ Int_t GetCurrentNodeId() const;
virtual TGeoNode *GetMother(Int_t up=1) const;
virtual TGeoNode *GetNode() const;
Int_t GetStackLevel() const {return fStackLevel;}
@@ -183,7 +185,7 @@ public:
virtual void LocalToMasterBomb(const Double_t *local, Double_t *master) const;
virtual void MasterToLocalBomb(const Double_t *master, Double_t *local) const;
virtual void PrintNode() const;
- virtual Int_t PushState(Bool_t ovlp, Double_t *point=0);
+ virtual Int_t PushState(Bool_t ovlp, Int_t startlevel=0, Double_t *point=0);
virtual Bool_t PopState(Double_t *point=0);
virtual Bool_t PopState(Int_t level, Double_t *point=0);
virtual void PopDummy(Int_t ipop=9999) {fStackLevel=(ipop>fStackLevel)?(fStackLevel-1):(ipop-1);}
@@ -209,10 +211,11 @@ private:
TGeoNode *fNode; // current node
TGeoHMatrix *fMatrix; // current matrix
TGeoHMatrix **fMatrixBranch; // current branch of global matrices
+ TGeoHMatrix **fMPB; // pre-built matrices
TGeoNode **fNodeBranch; // current branch of nodes
public:
TGeoCacheDummy();
- TGeoCacheDummy(TGeoNode *top);
+ TGeoCacheDummy(TGeoNode *top, Bool_t nodeid=kFALSE);
virtual ~TGeoCacheDummy();
virtual Bool_t CdDown(Int_t index, Bool_t make=kTRUE);
diff --git a/geom/inc/TGeoManager.h b/geom/inc/TGeoManager.h
index 04fe567258c..2718b850f1b 100644
--- a/geom/inc/TGeoManager.h
+++ b/geom/inc/TGeoManager.h
@@ -1,4 +1,4 @@
-// @(#)root/geom:$Name: $:$Id: TGeoManager.h,v 1.41 2003/12/10 15:31:23 brun Exp $
+// @(#)root/geom:$Name: $:$Id: TGeoManager.h,v 1.42 2004/01/19 13:40:50 brun Exp $
// Author: Andrei Gheata 25/10/01
/*************************************************************************
@@ -117,7 +117,7 @@ private :
Double_t *fDblBuffer; //! transient dbl buffer
//--- private methods
- void BuildCache(Bool_t dummy=kFALSE);
+ void BuildCache(Bool_t dummy=kFALSE, Bool_t nodeid=kFALSE);
void BuildIdArray();
TGeoNode *FindInCluster(Int_t *cluster, Int_t nc);
Int_t GetTouchedCluster(Int_t start, Double_t *point, Int_t *check_list,
@@ -310,6 +310,7 @@ public:
TVirtualGeoTrack *GetParentTrackOfId(Int_t id) const;
Int_t GetVirtualLevel();
Bool_t GotoSafeLevel();
+ Int_t GetSafeLevel() const;
Double_t GetSafeDistance() const {return fSafety;}
Double_t GetStep() const {return fStep;}
Bool_t IsAnimatingTracks() const {return fIsGeomReading;}
@@ -430,12 +431,12 @@ public:
void SelectTrackingMedia();
//--- stack manipulation
- Int_t PushPath() {return fCache->PushState(fCurrentOverlapping);}
+ Int_t PushPath(Int_t startlevel=0) {return fCache->PushState(fCurrentOverlapping, startlevel);}
Bool_t PopPath() {fCurrentOverlapping=fCache->PopState(); fCurrentNode=fCache->GetNode();
fLevel=fCache->GetLevel();return fCurrentOverlapping;}
Bool_t PopPath(Int_t index) {fCurrentOverlapping=fCache->PopState(index);
fCurrentNode=fCache->GetNode(); fLevel=fCache->GetLevel();return fCurrentOverlapping;}
- Int_t PushPoint() {return fCache->PushState(fCurrentOverlapping, fPoint);}
+ Int_t PushPoint(Int_t startlevel=0) {return fCache->PushState(fCurrentOverlapping, startlevel,fPoint);}
Bool_t PopPoint() {fCurrentOverlapping=fCache->PopState(fPoint); fCurrentNode=fCache->GetNode();
fLevel=fCache->GetLevel(); return fCurrentOverlapping;}
Bool_t PopPoint(Int_t index) {fCurrentOverlapping=fCache->PopState(index, fPoint); fCurrentNode=fCache->GetNode();
diff --git a/geom/inc/TGeoMatrix.h b/geom/inc/TGeoMatrix.h
index d782fd12c12..5f48f419d07 100644
--- a/geom/inc/TGeoMatrix.h
+++ b/geom/inc/TGeoMatrix.h
@@ -1,4 +1,4 @@
-// @(#)root/geom:$Name: $:$Id: TGeoMatrix.h,v 1.10 2003/09/25 14:50:40 brun Exp $
+// @(#)root/geom:$Name: $:$Id: TGeoMatrix.h,v 1.11 2004/01/18 12:31:54 brun Exp $
// Author: Andrei Gheata 25/10/01
/*************************************************************************
@@ -301,13 +301,14 @@ public :
TGeoIdentity(const char *name);
virtual ~TGeoIdentity() {}
- virtual void LocalToMaster(const Double_t *local, Double_t *master) const {memcpy(master, local, 3*sizeof(Double_t));}
+ virtual void LocalToMaster(const Double_t *local, Double_t *master) const
+ {memcpy(master, local, 3*sizeof(Double_t));}
virtual void LocalToMasterVect(const Double_t *local, Double_t *master) const
-{memcpy(master, local, 3*sizeof(Double_t));}
+ {memcpy(master, local, 3*sizeof(Double_t));}
virtual void MasterToLocal(const Double_t *master, Double_t *local) const
- {memcpy(local, master, 3*sizeof(Double_t));}
+ {memcpy(local, master, 3*sizeof(Double_t));}
virtual void MasterToLocalVect(const Double_t *master, Double_t *local) const
- {memcpy(local, master, 3*sizeof(Double_t));}
+ {memcpy(local, master, 3*sizeof(Double_t));}
virtual void LocalToMasterBomb(const Double_t *local, Double_t *master) const
{TGeoIdentity::LocalToMaster(local, master);}
virtual void MasterToLocalBomb(const Double_t *master, Double_t *local) const
diff --git a/geom/inc/TGeoVoxelFinder.h b/geom/inc/TGeoVoxelFinder.h
index 68c71f22e89..de4197edf79 100644
--- a/geom/inc/TGeoVoxelFinder.h
+++ b/geom/inc/TGeoVoxelFinder.h
@@ -1,4 +1,4 @@
-// @(#)root/geom:$Name: $:$Id: TGeoVoxelFinder.h,v 1.7 2003/02/18 15:37:36 brun Exp $
+// @(#)root/geom:$Name: $:$Id: TGeoVoxelFinder.h,v 1.8 2004/01/18 12:31:54 brun Exp $
// Author: Andrei Gheata 04/02/02
/*************************************************************************
@@ -23,6 +23,10 @@
class TGeoVoxelFinder : public TObject
{
+public:
+enum EVoxelsType {
+ kGeoInvalidVoxels = BIT(15)
+};
protected:
TGeoVolume *fVolume; // volume to which applies
@@ -72,6 +76,7 @@ public :
void DaughterToMother(Int_t id, Double_t *local, Double_t *master) const;
virtual Double_t Efficiency();
virtual Int_t *GetCheckList(Double_t *point, Int_t &nelem);
+ Int_t *GetCheckList(Int_t &nelem) const {nelem=fNcandidates; return fCheckList;}
Int_t *GetExtraX(Int_t islice, Bool_t left, Int_t &nextra) const;
Int_t *GetExtraY(Int_t islice, Bool_t left, Int_t &nextra) const;
Int_t *GetExtraZ(Int_t islice, Bool_t left, Int_t &nextra) const;
@@ -84,7 +89,9 @@ public :
Int_t *GetValidExtra(Int_t *list, Int_t &ncheck);
Int_t *GetValidExtra(Int_t n1, UChar_t *array1, Int_t *list, Int_t &ncheck);
Int_t *GetValidExtra(Int_t n1, UChar_t *array1, Int_t n2, UChar_t *array2, Int_t *list, Int_t &ncheck);
+ virtual Int_t *GetVoxelCandidates(Int_t i, Int_t j, Int_t k, Int_t &ncheck);
virtual void FindOverlaps(Int_t inode) const;
+ Bool_t IsInvalid() const {return TObject::TestBit(kGeoInvalidVoxels);}
Double_t *GetBoxes() const {return fBoxes;}
Bool_t IsSafeVoxel(Double_t *point, Int_t inode, Double_t minsafe) const;
virtual void Print(Option_t *option="") const;
@@ -102,6 +109,7 @@ public :
Bool_t IntersectAndStore(Int_t n1, UChar_t *array1, Int_t n2, UChar_t *array2);
Bool_t IntersectAndStore(Int_t n1, UChar_t *array1, Int_t n2, UChar_t *array2,
Int_t n3, UChar_t *array3);
+ void SetInvalid(Bool_t flag=kTRUE) {TObject::SetBit(kGeoInvalidVoxels, flag);}
virtual void SortAll(Option_t *option="");
void SortCrossedVoxels(Double_t *point, Double_t *dir);
// Bool_t Union(Int_t n1, Int_t *array1, Int_t n2, Int_t *array2,
@@ -143,4 +151,37 @@ public:
ClassDef(TGeoCylVoxels, 2) // cylindrical voxel class
};
+/*************************************************************************
+ * TGeoFullVoxels - Full voxeliztion applies to volumes with limited
+ * number of daughters and stores node information in each cell.
+ *
+ *************************************************************************/
+
+class TGeoFullVoxels : public TGeoVoxelFinder
+{
+private:
+ Int_t fNvoxels;
+ Int_t fNvx; // number of slices on X
+ Int_t fNvy; // number of slices on Y
+ Int_t fNvz; // number of slices on Z
+ UChar_t *fVox; //[fNvoxels] voxels storage array
+public:
+ TGeoFullVoxels();
+ TGeoFullVoxels(TGeoVolume *vol);
+ virtual ~TGeoFullVoxels();
+
+// virtual void BuildVoxelLimits();
+// virtual Double_t Efficiency();
+// virtual void FindOverlaps(Int_t inode) const;
+ virtual Int_t *GetCheckList(Double_t *point, Int_t &nelem);
+// virtual Bool_t GetNextIndices(Double_t *point, Double_t *dir);
+ UChar_t *GetVoxel(Int_t i, Int_t j, Int_t k) const {return &fVox[(i*fNvy+j)*fNvz+k];}
+ virtual Int_t *GetVoxelCandidates(Int_t i, Int_t j, Int_t k, Int_t &ncheck);
+// virtual Int_t *GetNextVoxel(Double_t *point, Double_t *dir, Int_t &ncheck);
+ virtual void Print(Option_t *option="") const;
+ virtual void Voxelize(Option_t *option);
+
+ ClassDef(TGeoFullVoxels, 1) // full voxels class
+};
+
#endif
diff --git a/geom/src/TGeoCache.cxx b/geom/src/TGeoCache.cxx
index 94c10bb1349..a99521ef6fa 100644
--- a/geom/src/TGeoCache.cxx
+++ b/geom/src/TGeoCache.cxx
@@ -1,4 +1,4 @@
-// @(#)root/geom:$Name: $:$Id: TGeoCache.cxx,v 1.23 2004/01/19 13:43:15 brun Exp $
+// @(#)root/geom:$Name: $:$Id: TGeoCache.cxx,v 1.24 2004/01/19 13:44:14 brun Exp $
// Author: Andrei Gheata 18/03/02
/*************************************************************************
@@ -68,11 +68,45 @@ TGeoNodeCache::TGeoNodeCache()
fMatrixPool = 0;
fNodeIdArray = 0;
fIndex = 0;
- BuildIdArray();
+// BuildIdArray();
}
//_____________________________________________________________________________
-TGeoNodeCache::TGeoNodeCache(Int_t /*size*/)
+TGeoNodeCache::TGeoNodeCache(Bool_t nodeid)
+{
+// dummy constructor
+ fGeoCacheMaxDaughters = 128;
+ fGeoCacheMaxSize = 1000000;
+ fGeoCacheStackSize = 1000;
+ fGeoCacheDefaultLevel = 4;
+ fGeoCacheMaxLevels = 30;
+ fGeoCacheObjArrayInd = 0xFF;
+ fGeoCacheUsageRatio = 0.01;
+ fSize = 0;
+ fNused = 0;
+ fLevel = 0;
+ fStackLevel = 0;
+ fStack = 0;
+ fDefaultLevel= 0;
+ fCache = 0;
+ fPath = "";
+ fTopNode = 0;
+ fCurrentNode = 0;
+ fCurrentCache = 0;
+ fCurrentIndex = 0;
+ fCurrentID = 0;
+ fBranch = 0;
+ fMatrices = 0;
+ fGlobalMatrix= 0;
+ fMatrixPool = 0;
+ fNodeIdArray = 0;
+ fIndex = 0;
+ if (nodeid) BuildIdArray();
+ else printf("--- node ID tracking disabled\n");
+}
+
+//_____________________________________________________________________________
+TGeoNodeCache::TGeoNodeCache(Int_t /*size*/, Bool_t nodeid)
{
// constructor
fGeoCacheMaxDaughters = 128;
@@ -120,7 +154,8 @@ TGeoNodeCache::TGeoNodeCache(Int_t /*size*/)
fCurrentID = 0;
fNodeIdArray = 0;
fIndex = 0;
- BuildIdArray();
+ if (nodeid) BuildIdArray();
+ else printf("--- node ID tracking disabled\n");
CdTop();
}
@@ -148,8 +183,8 @@ void TGeoNodeCache::BuildIdArray()
// Builds node id array.
Int_t nnodes = gGeoManager->GetNNodes();
//if (nnodes>3E7) return;
- if (nnodes>0) return; //do not use the node cache by default
if (fNodeIdArray) delete [] fNodeIdArray;
+ printf("--- node ID tracking enabled, size=%i Bytes\n", (2*nnodes+1)*sizeof(Int_t));
fNodeIdArray = new Int_t[2*nnodes+1];
fNodeIdArray[0] = 0;
Int_t ifree = 1;
@@ -158,6 +193,14 @@ void TGeoNodeCache::BuildIdArray()
fIdBranch[0] = 0;
}
+//_____________________________________________________________________________
+Int_t TGeoNodeCache::GetCurrentNodeId() const
+{
+// Returns a fixed ID for current physical node
+ if (fNodeIdArray) return fNodeIdArray[fIndex];
+ return GetNodeId();
+}
+
//_____________________________________________________________________________
void TGeoNodeCache::Compact()
{
@@ -469,13 +512,13 @@ void TGeoNodeCache::PrintNode() const
}
//_____________________________________________________________________________
-Int_t TGeoNodeCache::PushState(Bool_t ovlp, Double_t *point)
+Int_t TGeoNodeCache::PushState(Bool_t ovlp, Int_t startlevel, Double_t *point)
{
if (fStackLevel>=fGeoCacheStackSize) {
printf("ERROR TGeoNodeCach::PushSate() : stack of states full\n");
return 0;
}
- ((TGeoCacheState*)fStack->At(fStackLevel))->SetState(fLevel,ovlp,point);
+ ((TGeoCacheState*)fStack->At(fStackLevel))->SetState(fLevel,startlevel,ovlp,point);
return ++fStackLevel;
}
@@ -549,19 +592,24 @@ TGeoCacheDummy::TGeoCacheDummy()
fNode = 0;
fNodeBranch = 0;
fMatrixBranch = 0;
+ fMPB = 0;
}
//_____________________________________________________________________________
-TGeoCacheDummy::TGeoCacheDummy(TGeoNode *top)
+TGeoCacheDummy::TGeoCacheDummy(TGeoNode *top, Bool_t nodeid)
+ :TGeoNodeCache(nodeid)
{
fTop = top;
fNode = top;
fNodeBranch = new TGeoNode *[fGeoCacheMaxLevels];
fNodeBranch[0] = top;
fMatrixBranch = new TGeoHMatrix *[fGeoCacheMaxLevels];
- for (Int_t i=0; i<fGeoCacheMaxLevels; i++)
- fMatrixBranch[i] = new TGeoHMatrix("global");
- fMatrix = fMatrixBranch[0];
+ fMPB = new TGeoHMatrix *[fGeoCacheMaxLevels];
+ for (Int_t i=0; i<fGeoCacheMaxLevels; i++) {
+ fMPB[i] = new TGeoHMatrix("global");
+ fMatrixBranch[i] = 0;
+ }
+ fMatrix = fMatrixBranch[0] = fMPB[0];
fStack = new TObjArray(fGeoCacheStackSize);
for (Int_t ist=0; ist<fGeoCacheStackSize; ist++)
fStack->Add(new TGeoCacheStateDummy(100)); // !obsolete 100
@@ -572,10 +620,10 @@ TGeoCacheDummy::TGeoCacheDummy(TGeoNode *top)
TGeoCacheDummy::~TGeoCacheDummy()
{
if (fNodeBranch) delete [] fNodeBranch;
- if (fMatrixBranch) {
+ if (fMPB) {
for (Int_t i=0; i<fGeoCacheMaxLevels; i++)
- delete fMatrixBranch[i];
- delete [] fMatrixBranch;
+ delete fMPB[i];
+ delete [] fMPB;
}
}
@@ -584,18 +632,20 @@ Bool_t TGeoCacheDummy::CdDown(Int_t index, Bool_t /*make*/)
{
TGeoNode *newnode = fNode->GetDaughter(index);
if (!newnode) return kFALSE;
- TGeoHMatrix *newmat = fMatrixBranch[fLevel+1];
- TGeoMatrix *local = newnode->GetMatrix();
- *newmat = fMatrix;
- newmat->Multiply(local);
fLevel++;
if (fNodeIdArray) {
fIndex = fNodeIdArray[fIndex+index+1];
fIdBranch[fLevel] = fIndex;
}
- fMatrix = newmat;
fNode = newnode;
fNodeBranch[fLevel] = fNode;
+ TGeoMatrix *local = newnode->GetMatrix();
+ TGeoHMatrix *newmat = fMPB[fLevel+1];
+ if (!local->IsIdentity()) {
+ *newmat = fMatrix;
+ newmat->Multiply(local);
+ fMatrix = newmat;
+ }
fMatrixBranch[fLevel] = fMatrix;
return kTRUE;
}
@@ -1457,17 +1507,20 @@ TGeoCacheState::~TGeoCacheState()
}
//_____________________________________________________________________________
-void TGeoCacheState::SetState(Int_t level, Bool_t ovlp, Double_t *point)
+void TGeoCacheState::SetState(Int_t level, Int_t startlevel, Bool_t ovlp, Double_t *point)
{
fLevel = level;
+ fStart = startlevel;
if (gGeoManager->IsOutside()) {
fLevel = -1;
return;
}
TGeoNodeCache *cache = gGeoManager->GetCache();
- if (cache->HasIdArray()) memcpy(fIdBranch, cache->GetIdBranch(), (level+1)*sizeof(Int_t));
- memcpy(fBranch, (Int_t*)cache->GetBranch(), (level+1)*sizeof(Int_t));
- memcpy(fMatrices, (Int_t*)cache->GetMatrices(), (level+1)*sizeof(Int_t));
+ Int_t *branch = (Int_t*)cache->GetBranch();
+ Int_t *matrices = (Int_t*)cache->GetMatrices();
+ if (cache->HasIdArray()) memcpy(fIdBranch, cache->GetIdBranch()+fStart, (level+1-fStart)*sizeof(Int_t));
+ memcpy(fBranch, branch+fStart, (level+1-fStart)*sizeof(Int_t));
+ memcpy(fMatrices, matrices+fStart, (level+1-fStart)*sizeof(Int_t));
fOverlapping = ovlp;
if (point) memcpy(fPoint, point, 3*sizeof(Double_t));
}
@@ -1481,9 +1534,11 @@ Bool_t TGeoCacheState::GetState(Int_t &level, Double_t *point) const
return kFALSE;
}
TGeoNodeCache *cache = gGeoManager->GetCache();
- if (cache->HasIdArray()) cache->FillIdBranch(fIdBranch);
- memcpy((Int_t*)cache->GetBranch(), fBranch, (level+1)*sizeof(Int_t));
- memcpy((Int_t*)cache->GetMatrices(), fMatrices, (level+1)*sizeof(Int_t));
+ if (cache->HasIdArray()) cache->FillIdBranch(fIdBranch, fStart);
+ Int_t *branch = (Int_t*)cache->GetBranch();
+ Int_t *matrices = (Int_t*)cache->GetMatrices();
+ memcpy(branch+fStart, fBranch, (level+1-fStart)*sizeof(Int_t));
+ memcpy(matrices+fStart, fMatrices, (level+1-fStart)*sizeof(Int_t));
if (point) memcpy(point, fPoint, 3*sizeof(Double_t));
return fOverlapping;
}
@@ -1531,17 +1586,18 @@ TGeoCacheStateDummy::~TGeoCacheStateDummy()
}
//_____________________________________________________________________________
-void TGeoCacheStateDummy::SetState(Int_t level, Bool_t ovlp, Double_t *point)
+void TGeoCacheStateDummy::SetState(Int_t level, Int_t startlevel, Bool_t ovlp, Double_t *point)
{
fLevel = level;
+ fStart = startlevel;
TGeoNodeCache *cache = gGeoManager->GetCache();
- if (cache->HasIdArray()) memcpy(fIdBranch, cache->GetIdBranch(), (level+1)*sizeof(Int_t));
+ if (cache->HasIdArray()) memcpy(fIdBranch, cache->GetIdBranch()+fStart, (level+1-fStart)*sizeof(Int_t));
TGeoNode **node_branch = (TGeoNode **) cache->GetBranch();
TGeoHMatrix **mat_branch = (TGeoHMatrix **) cache->GetMatrices();
- memcpy(fNodeBranch, node_branch, (level+1)*sizeof(TGeoNode *));
- for (Int_t i=0; i<level+1; i++)
- *fMatrixBranch[i] = mat_branch[i];
+ memcpy(fNodeBranch, node_branch+fStart, (level+1-fStart)*sizeof(TGeoNode *));
+ for (Int_t i=0; i<level+1-fStart; i++)
+ *fMatrixBranch[i] = mat_branch[i+fStart];
fOverlapping = ovlp;
if (point) memcpy(fPoint, point, 3*sizeof(Double_t));
}
@@ -1551,13 +1607,13 @@ Bool_t TGeoCacheStateDummy::GetState(Int_t &level, Double_t *point) const
{
level = fLevel;
TGeoNodeCache *cache = gGeoManager->GetCache();
- if (cache->HasIdArray()) cache->FillIdBranch(fIdBranch);
+ if (cache->HasIdArray()) cache->FillIdBranch(fIdBranch, fStart);
TGeoNode **node_branch = (TGeoNode **) cache->GetBranch();
TGeoHMatrix **mat_branch = (TGeoHMatrix **) cache->GetMatrices();
- memcpy(node_branch, fNodeBranch, (level+1)*sizeof(TGeoNode *));
- for (Int_t i=0; i<level+1; i++)
- *mat_branch[i] = fMatrixBranch[i];
+ memcpy(node_branch+fStart, fNodeBranch, (level+1-fStart)*sizeof(TGeoNode *));
+ for (Int_t i=0; i<level+1-fStart; i++)
+ *mat_branch[i+fStart] = fMatrixBranch[i];
if (point) memcpy(point, fPoint, 3*sizeof(Double_t));
return fOverlapping;
}
diff --git a/geom/src/TGeoManager.cxx b/geom/src/TGeoManager.cxx
index 03bd9cd44ab..10614e3a67a 100644
--- a/geom/src/TGeoManager.cxx
+++ b/geom/src/TGeoManager.cxx
@@ -1,4 +1,4 @@
-// @(#)root/geom:$Name: $:$Id: TGeoManager.cxx,v 1.69 2004/01/18 12:31:54 brun Exp $
+// @(#)root/geom:$Name: $:$Id: TGeoManager.cxx,v 1.70 2004/01/19 13:40:51 brun Exp $
// Author: Andrei Gheata 25/10/01
/*************************************************************************
@@ -714,16 +714,16 @@ void TGeoManager::UnbombTranslation(const Double_t *tr, Double_t *bombtr)
}
//_____________________________________________________________________________
-void TGeoManager::BuildCache(Bool_t dummy)
+void TGeoManager::BuildCache(Bool_t dummy, Bool_t nodeid)
{
// Builds the cache for physical nodes and global matrices.
if (!fCache) {
if (fNNodes>5000000 || dummy) // temporary - works without
// build dummy cache
- fCache = new TGeoCacheDummy(fTopNode);
+ fCache = new TGeoCacheDummy(fTopNode, nodeid);
else
// build real cache
- fCache = new TGeoNodeCache(0);
+ fCache = new TGeoNodeCache(0,nodeid);
}
}
@@ -1333,6 +1333,7 @@ void TGeoManager::CloseGeometry(Option_t *option)
TString opt(option);
opt.ToLower();
Bool_t dummy = opt.Contains("d");
+ Bool_t nodeid = opt.Contains("i");
if (fIsGeomReading) {
printf("### Geometry loaded from file...\n");
gGeoIdentity=(TGeoIdentity *)fMatrices->At(0);
@@ -1344,11 +1345,11 @@ void TGeoManager::CloseGeometry(Option_t *option)
SetTopVolume(fMasterVolume);
if (fStreamVoxels) printf("### Voxelization retrieved from file\n");
Voxelize("ALL");
- if (!fCache) BuildCache(dummy);
+ if (!fCache) BuildCache(dummy,nodeid);
} else {
Warning("CloseGeometry", "top node was streamed!");
Voxelize("ALL");
- if (!fCache) BuildCache(dummy);
+ if (!fCache) BuildCache(dummy,nodeid);
}
// BuildIdArray();
printf("### %i nodes/ %i volume UID's in %s\n", fNNodes, fUniqueVolumes->GetEntriesFast(), GetTitle());
@@ -1364,7 +1365,7 @@ void TGeoManager::CloseGeometry(Option_t *option)
// BuildIdArray();
Voxelize("ALL");
printf("Building caches for nodes and matrices...\n");
- BuildCache(dummy);
+ BuildCache(dummy,nodeid);
printf("### %i nodes/ %i volume UID's in %s\n", fNNodes, fUniqueVolumes->GetEntriesFast(), GetTitle());
printf("----------------modeler ready----------------\n");
}
@@ -1766,6 +1767,23 @@ Bool_t TGeoManager::GotoSafeLevel()
while (fCurrentOverlapping && fLevel) CdUp();
return kTRUE;
}
+
+//_____________________________________________________________________________
+Int_t TGeoManager::GetSafeLevel() const
+{
+// Go upwards the tree until a non-overlaping node
+ Bool_t overlapping = fCurrentOverlapping;
+ if (!overlapping) return fLevel;
+ Int_t level = fLevel;
+ TGeoNode *node;
+ while (overlapping && level) {
+ level--;
+ node = GetMother(fLevel-level);
+ if (!node->IsOffset()) overlapping = node->IsOverlapping();
+ }
+ return level;
+}
+
//_____________________________________________________________________________
TGeoNode *TGeoManager::FindInCluster(Int_t *cluster, Int_t nc)
{
@@ -2018,7 +2036,8 @@ void TGeoManager::SafetyOverlaps()
TGeoVolume *vol;
Int_t novlp, io;
Int_t *ovlp;
- PushPath();
+ Int_t safelevel = GetSafeLevel();
+ PushPath(safelevel+1);
while (fCurrentOverlapping) {
ovlp = fCurrentNode->GetOverlaps(novlp);
CdUp();
@@ -2571,9 +2590,10 @@ TGeoNode *TGeoManager::FindNextBoundary(Double_t stepmax, const char *path)
Double_t mothpt[3];
Double_t vecpt[3];
Double_t dpt[3], dvec[3];
- PushPath();
Int_t novlps;
- while (fCurrentOverlapping) {
+ Int_t safelevel = GetSafeLevel();
+ PushPath(safelevel+1);
+ while (fCurrentOverlapping) {
Int_t *ovlps = fCurrentNode->GetOverlaps(novlps);
CdUp();
mother = fCurrentNode->GetVolume();
@@ -2885,13 +2905,10 @@ Bool_t TGeoManager::IsSameLocation(Double_t x, Double_t y, Double_t z, Bool_t ch
// Checks if point (x,y,z) is still in the current node.
// check if this is an overlapping node
if (fCurrentOverlapping) {
- TGeoNode *current = fCurrentNode;
Int_t cid = GetCurrentNodeId();
if (!change) PushPoint();
gGeoManager->FindNode(x,y,z);
- Bool_t same = (cid>=0)?kTRUE:kFALSE;
- if (same) same = (cid==GetCurrentNodeId())?kTRUE:kFALSE;
- else same = (current==fCurrentNode)?kTRUE:kFALSE;
+ Bool_t same = (cid==GetCurrentNodeId())?kTRUE:kFALSE;
if (!change) PopPoint();
return same;
}
@@ -3561,9 +3578,10 @@ void TGeoManager::SetTopVolume(TGeoVolume *vol)
fLevel = 0;
if (fCache) {
Bool_t dummy=fCache->IsDummy();
+ Bool_t nodeid = fCache->HasIdArray();
delete fCache;
fCache = 0;
- BuildCache(dummy);
+ BuildCache(dummy,nodeid);
}
printf("Top volume is %s. Master volume is %s\n", fTopVolume->GetName(),
fMasterVolume->GetName());
diff --git a/geom/src/TGeoTube.cxx b/geom/src/TGeoTube.cxx
index c99ff1c7aef..78f9520b4e7 100644
--- a/geom/src/TGeoTube.cxx
+++ b/geom/src/TGeoTube.cxx
@@ -1,4 +1,4 @@
-// @(#)root/geom:$Name: $:$Id: TGeoTube.cxx,v 1.29 2003/12/10 15:31:23 brun Exp $
+// @(#)root/geom:$Name: $:$Id: TGeoTube.cxx,v 1.30 2003/12/11 10:34:33 brun Exp $
// Author: Andrei Gheata 24/10/01
// TGeoTube::Contains() and DistToOut/In() implemented by Mihaela Gheata
@@ -303,7 +303,7 @@ Double_t TGeoTube::DistToOut(Double_t *point, Double_t *dir, Int_t iact, Double_
Double_t b,d;
Double_t sr;
// inner cylinder
- if (fRmin>1E-10) {
+ if (fRmin>1E-10 && rdotn<0) {
DistToTube(rsq,nsq,rdotn,fRmin,b,d);
if (d>0) {
sr=-b-d;
diff --git a/geom/src/TGeoVolume.cxx b/geom/src/TGeoVolume.cxx
index 0e3ea82403e..ff09ed631e3 100644
--- a/geom/src/TGeoVolume.cxx
+++ b/geom/src/TGeoVolume.cxx
@@ -1,4 +1,4 @@
-// @(#)root/geom:$Name: $:$Id: TGeoVolume.cxx,v 1.39 2003/10/20 08:46:33 brun Exp $
+// @(#)root/geom:$Name: $:$Id: TGeoVolume.cxx,v 1.40 2003/11/11 15:44:28 brun Exp $
// Author: Andrei Gheata 30/05/02
// Divide(), CheckOverlaps() implemented by Mihaela Gheata
@@ -1346,15 +1346,20 @@ void TGeoVolume::Voxelize(Option_t *option)
}
}
// find optimal voxelization
- Bool_t cyltype = GetOptimalVoxels();
- if (cyltype) {
-// fVoxels = new TGeoCylVoxels(this);
+// Bool_t cyltype = GetOptimalVoxels();
+// if (nd < 8) {
+// fVoxels = new TGeoFullVoxels(this);
+// if (fVoxels) printf("full voxels for %s\n", GetName());
+// } else
fVoxels = new TGeoVoxelFinder(this);
-// printf("%s cyl. voxels\n", GetName());
- } else {
- fVoxels = new TGeoVoxelFinder(this);
- }
+
fVoxels->Voxelize(option);
+ if (fVoxels) {
+ if (fVoxels->IsInvalid()) {
+ delete fVoxels;
+ fVoxels = 0;
+ }
+ }
// if (fVoxels) fVoxels->Print();
}
diff --git a/geom/src/TGeoVoxelFinder.cxx b/geom/src/TGeoVoxelFinder.cxx
index d72836cab6e..40f5a40937f 100644
--- a/geom/src/TGeoVoxelFinder.cxx
+++ b/geom/src/TGeoVoxelFinder.cxx
@@ -1,4 +1,4 @@
-// @(#)root/geom:$Name: $:$Id: TGeoVoxelFinder.cxx,v 1.20 2003/12/11 10:34:33 brun Exp $
+// @(#)root/geom:$Name: $:$Id: TGeoVoxelFinder.cxx,v 1.21 2004/01/18 12:31:55 brun Exp $
// Author: Andrei Gheata 04/02/02
/*************************************************************************
@@ -75,6 +75,7 @@ TGeoVoxelFinder::TGeoVoxelFinder()
fCheckList = 0;
fNcandidates = 0;
fCurrentVoxel = 0;
+ SetInvalid(kFALSE);
}
//-----------------------------------------------------------------------------
TGeoVoxelFinder::TGeoVoxelFinder(TGeoVolume *vol)
@@ -1129,19 +1130,27 @@ void TGeoVoxelFinder::SortCrossedVoxels(Double_t *point, Double_t *dir)
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;
- Bool_t one_daughter = kFALSE;
- Int_t nslices = 0;
+// 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;
- one_daughter = kTRUE;
+ return fCheckList;
}
+ Int_t nslices = 0;
UChar_t *slice1 = 0;
UChar_t *slice2 = 0;
UChar_t *slice3 = 0;
- Int_t nd[3];
- memset(&nd[0], 0, 3*sizeof(Int_t));
+ Int_t nd[3] = {0,0,0};
Int_t im;
if (fPriority[0]) {
im = TMath::BinarySearch(fIbx, fXb, point[0]);
@@ -1190,7 +1199,6 @@ Int_t *TGeoVoxelFinder::GetCheckList(Double_t *point, Int_t &nelem)
}
}
}
- if (one_daughter) return fCheckList;
nelem = 0;
// Int_t i = 0;
Bool_t intersect = kFALSE;
@@ -1210,6 +1218,69 @@ Int_t *TGeoVoxelFinder::GetCheckList(Double_t *point, Int_t &nelem)
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)
{
@@ -1486,7 +1557,6 @@ void TGeoVoxelFinder::SortAll(Option_t *)
{
// order bounding boxes along x, y, z
Int_t nd = fVolume->GetNdaughters();
- if (!nd) return;
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
@@ -1545,13 +1615,14 @@ void TGeoVoxelFinder::SortAll(Option_t *)
}
// now find priority
if (ib < 2) {
- Error("SortAll", "less than 2 boundaries on X for %s!", fVolume->GetName());
+ 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) {
@@ -1663,13 +1734,14 @@ void TGeoVoxelFinder::SortAll(Option_t *)
}
// now find priority on Y
if (ib < 2) {
- Error("SortAll", "less than 2 boundaries on Y for %s!", fVolume->GetName());
+ 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) {
@@ -1783,13 +1855,14 @@ void TGeoVoxelFinder::SortAll(Option_t *)
}
// now find priority on Z
if (ib < 2) {
- Error("SortAll", "less than 2 boundaries on Z for %s!", fVolume->GetName());
+ 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) {
@@ -1819,7 +1892,6 @@ void TGeoVoxelFinder::SortAll(Option_t *)
fZb = new Double_t[ib];
memcpy(fZb, &temp[0], ib*sizeof(Double_t));
fIbz = ib;
-
}
@@ -2522,6 +2594,7 @@ void TGeoCylVoxels::SortAll(Option_t *)
if ((rmin>=rmax) || (pmin>=pmax) || (zmin>=zmax)) {
Error("SortAll", "wrong bounding cylinder");
printf("### volume was : %s\n", fVolume->GetName());
+ SetInvalid();
return;
}
Int_t id;
@@ -2572,6 +2645,7 @@ void TGeoCylVoxels::SortAll(Option_t *)
if (ib < 2) {
Error("SortAll", "less than 2 boundaries on R !");
printf("### volume was : %s\n", fVolume->GetName());
+ SetInvalid();
return;
}
if (ib == 2) {
@@ -2654,6 +2728,7 @@ void TGeoCylVoxels::SortAll(Option_t *)
if (ib < 2) {
Error("SortAll", "less than 2 boundaries on Phi !");
printf("### volume was : %s\n", fVolume->GetName());
+ SetInvalid();
return;
}
if (ib == 2) {
@@ -2726,6 +2801,7 @@ void TGeoCylVoxels::SortAll(Option_t *)
if (ib < 2) {
Error("SortAll", "less than 2 boundaries on Z !");
printf("### volume was : %s\n", fVolume->GetName());
+ SetInvalid();
return;
}
if (ib == 2) {
@@ -2810,8 +2886,184 @@ void TGeoCylVoxels::Voxelize(Option_t *)
{
//--- Voxelize fVolume.
// printf("Voxelizing %s\n", fVolume->GetName());
- Int_t nd = fVolume->GetNdaughters();
- if (!nd) return;
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
+{
+}
+
--
GitLab