| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 #include "OctreeNode.hh"
28
29 #include "G4VPhysicalVolume.hh"
30
31 #include <cassert>
32 #include <utility>
33
34 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
35
36 OctreeNode::OctreeNode(const G4ThreeVector& position, const G4ThreeVector& halfLengths,
37 G4int maxContents, OctreeNode* parent)
38 : fPosition(position),
39 fHalfLengths(halfLengths),
40 fMaxContents(maxContents),
41 fParent(parent),
42 fChildren(std::array<OctreeNode*, 8>())
43 {
44 fHalfLengthsMag = fHalfLengths.mag();
45 }
46
47 OctreeNode::~OctreeNode()
48 {
49 // Delete children
50 for (auto& it : fChildren) {
51 delete it;
52 }
53 }
54
55 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
56
57 const OctreeNode* OctreeNode::GetChildFromPosition(const G4ThreeVector& pos) const
58 {
59 G4ThreeVector localPosition = pos - fPosition;
60
61 // Get the right child
62 // The scheme is defined by quadrant as follows:
63 // Zero is treated as positive
64 // X Y Z Index | X Y Z Index
65 // + + + 0 | - + + 4
66 // + + - 1 | - + - 5
67 // + - + 2 | - - + 6
68 // + - - 3 | - - - 7
69 int bit2 = (localPosition.getX() < 0) ? 4 : 0;
70 int bit1 = (localPosition.getY() < 0) ? 2 : 0;
71 int bit0 = (localPosition.getZ() < 0) ? 1 : 0;
72 return fChildren[bit0 + bit1 + bit2];
73 }
74
75 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
76
77 // Non-const version of above function. Ugly but useful in this case
78 OctreeNode* OctreeNode::GetChildFromPosition(const G4ThreeVector& pos)
79 {
80 G4ThreeVector localPosition = pos - fPosition;
81
82 // Get the right child
83 // The scheme is defined by quadrant as follows:
84 // Zero is treated as positive
85 // X Y Z Index | X Y Z Index
86 // + + + 0 | - + + 4
87 // + + - 1 | - + - 5
88 // + - + 2 | - - + 6
89 // + - - 3 | - - - 7
90 int bit2 = (localPosition.getX() < 0) ? 4 : 0;
91 int bit1 = (localPosition.getY() < 0) ? 2 : 0;
92 int bit0 = (localPosition.getZ() < 0) ? 1 : 0;
93 int idx = bit0 + bit1 + bit2;
94
95 assert(idx < (int)fChildren.size());
96 return fChildren[idx];
97 }
98
99 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
100
101 void OctreeNode::AddPhysicalVolume(G4VPhysicalVolume* pv)
102 {
103 // Consider adding test for the validity of the PV position
104 if (this->HasChildren()) {
105 OctreeNode* child = this->GetChildFromPosition(pv->GetTranslation());
106 child->AddPhysicalVolume(pv);
107 }
108 else {
109 // if there are fMaxContents elements in the bin, we need to split
110 // the bin before adding a new element.
111 if ((G4int)fContents.size() == fMaxContents) {
112 this->Split();
113 this->AddPhysicalVolume(pv);
114 }
115 else {
116 fContents.push_back(pv);
117 }
118 }
119 }
120
121 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
122
123 void OctreeNode::Split()
124 {
125 G4ThreeVector hl = 0.5 * fHalfLengths;
126 G4double xp = fPosition.getX();
127 G4double yp = fPosition.getY();
128 G4double zp = fPosition.getZ();
129 G4double xl = hl.getX();
130 G4double yl = hl.getY();
131 G4double zl = hl.getZ();
132 G4ThreeVector newpos;
133
134 // Construct children
135 newpos = G4ThreeVector(xp + xl, yp + yl, zp + zl);
136 fChildren[0] = new OctreeNode(newpos, hl, fMaxContents, this);
137 newpos = G4ThreeVector(xp + xl, yp + yl, zp - zl);
138 fChildren[1] = new OctreeNode(newpos, hl, fMaxContents, this);
139 newpos = G4ThreeVector(xp + xl, yp - yl, zp + zl);
140 fChildren[2] = new OctreeNode(newpos, hl, fMaxContents, this);
141 newpos = G4ThreeVector(xp + xl, yp - yl, zp - zl);
142 fChildren[3] = new OctreeNode(newpos, hl, fMaxContents, this);
143 newpos = G4ThreeVector(xp - xl, yp + yl, zp + zl);
144 fChildren[4] = new OctreeNode(newpos, hl, fMaxContents, this);
145 newpos = G4ThreeVector(xp - xl, yp + yl, zp - zl);
146 fChildren[5] = new OctreeNode(newpos, hl, fMaxContents, this);
147 newpos = G4ThreeVector(xp - xl, yp - yl, zp + zl);
148 fChildren[6] = new OctreeNode(newpos, hl, fMaxContents, this);
149 newpos = G4ThreeVector(xp - xl, yp - yl, zp - zl);
150 fChildren[7] = new OctreeNode(newpos, hl, fMaxContents, this);
151
152 // Distribute contents to children
153
154 for (int i = fContents.size() - 1; i >= 0; --i) {
155 G4VPhysicalVolume* pv = fContents[i];
156 OctreeNode* child = this->GetChildFromPosition(pv->GetTranslation());
157 assert(child != this);
158 child->AddPhysicalVolume(pv);
159 }
160 fContents.clear();
161 }
162
163 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
164
165 std::vector<G4VPhysicalVolume*> OctreeNode::GetContents() const
166 {
167 if (this->HasChildren()) // if has sub-nodes
168 {
169 std::vector<G4VPhysicalVolume*> vec;
170 std::vector<G4VPhysicalVolume*> childCont;
171 for (auto it = fChildren.begin(); it != fChildren.end(); ++it) {
172 childCont = (*it)->GetContents();
173 for (auto jt = childCont.begin(); jt != childCont.end(); ++jt) {
174 vec.push_back(*jt);
175 }
176 }
177 return vec;
178 }
179 else // if no sub-nodes
180 {
181 return fContents;
182 }
183 }
184
185 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
186
187 // Search for Octree nodes in a volume
188 const std::vector<G4VPhysicalVolume*> OctreeNode::SearchOctree(const G4ThreeVector& pos,
189 G4double _rad) const
190 {
191 // Need to search based on absolute position of each volume rather than
192 // relative position, and maintain a list of candidate nodes and
193 // final nodes within the radius
194
195 std::vector<const OctreeNode*> nodes;
196 std::vector<const OctreeNode*> candidates;
197 nodes.reserve(512);
198 candidates.reserve(512);
199
200 if (this->HasChildren()) {
201 for (auto it = this->GetChildren().begin(); it != this->GetChildren().end(); it++) {
202 candidates.push_back(*it);
203 }
204 }
205 else {
206 candidates.push_back(this);
207 }
208
209 const OctreeNode* aNode;
210 G4double d;
211 while (!candidates.empty()) {
212 aNode = candidates.back();
213 candidates.pop_back();
214 d = (aNode->GetPosition() - pos).mag() - aNode->GetHalfLengthsMag();
215 // if node within circle
216 if (d < _rad) {
217 if (aNode->HasChildren()) {
218 for (auto it = aNode->GetChildren().begin(); it != aNode->GetChildren().end(); ++it) {
219 candidates.push_back(*it);
220 }
221 }
222 else {
223 nodes.push_back(aNode);
224 }
225 }
226 }
227
228 // G4cout << "Found " << nodes.size() << " nodes" << G4endl;
229
230 // Get the physical volumes
231 G4double r2 = _rad * _rad;
232 std::vector<G4VPhysicalVolume*> pvols;
233 for (auto it = nodes.begin(); it != nodes.end(); ++it) {
234 std::vector<G4VPhysicalVolume*> conts = (*it)->GetContents();
235 for (auto jt = conts.begin(); jt != conts.end(); ++jt) {
236 if ((pos - (*jt)->GetTranslation()).mag2() < r2) {
237 pvols.push_back((*jt));
238 }
239 }
240 }
241 // G4cout << "Found " << pvols.size() << " candidates" << G4endl;
242 return pvols;
243 }
244
245 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
246
247 // Search for Octree nodes in a volume
248 void OctreeNode::SearchOctree(const G4ThreeVector& pos, std::vector<G4VPhysicalVolume*>& out,
249 G4double _rad) const
250 {
251 // Need to search based on absolute position of each volume rather than
252 // relative position, and maintain a list of candidate nodes and
253 // final nodes within the radius
254
255 std::vector<const OctreeNode*> nodes;
256 std::vector<const OctreeNode*> candidates;
257 nodes.reserve(512);
258 candidates.reserve(512);
259
260 if (this->HasChildren()) {
261 for (auto it = this->GetChildren().begin(); it != this->GetChildren().end(); ++it) {
262 candidates.push_back(*it);
263 }
264 }
265 else {
266 candidates.push_back(this);
267 }
268
269 const OctreeNode* aNode;
270 G4double d;
271 while (!candidates.empty()) {
272 aNode = candidates.back();
273 candidates.pop_back();
274 d = (aNode->GetPosition() - pos).mag() - aNode->GetHalfLengthsMag();
275 // if node within circle
276 if (d < _rad) {
277 if (aNode->HasChildren()) {
278 for (auto it = aNode->GetChildren().begin(); it != aNode->GetChildren().end(); it++) {
279 candidates.push_back(*it);
280 }
281 }
282 else {
283 nodes.push_back(aNode);
284 }
285 }
286 }
287
288 // Get the physical volumes
289 G4double r2 = _rad * _rad;
290 for (auto it = nodes.begin(); it != nodes.end(); ++it) {
291 std::vector<G4VPhysicalVolume*> conts = (*it)->GetContents();
292 for (auto jt = conts.begin(); jt != conts.end(); ++jt) {
293 if ((pos - (*jt)->GetTranslation()).mag2() < r2) {
294 out.push_back((*jt));
295 }
296 }
297 }
298 }
299
300 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
301
302 const std::vector<G4VPhysicalVolume*> OctreeNode::SearchOctree(const G4ThreeVector& pos) const
303 {
304 const OctreeNode* child = GetChildFromPosition(pos);
305 return child->GetContents();
306 }
307
308 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
309
310 G4int OctreeNode::GetNumberOfTerminalNodes()
311 {
312 if (!this->HasChildren()) {
313 return 1;
314 }
315 // sum the number of nodes of each child
316 G4int n = 0;
317 for (auto it = this->GetChildren().begin(); it != this->GetChildren().end(); ++it) {
318 n = n + (*it)->GetNumberOfTerminalNodes();
319 }
320
321 return n;
322 }
323 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
324