| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27 #include "OctreeNode.hh"
28
29 #include "G4VPhysicalVolume.hh"
30
31 #include <cassert>
32 #include <utility>
33
34 //....oooOO0OOooo........oooOO0OOooo........oo
35
36 OctreeNode::OctreeNode(const G4ThreeVector& po
37 G4int maxContents, Octr
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........oo
56
57 const OctreeNode* OctreeNode::GetChildFromPosi
58 {
59 G4ThreeVector localPosition = pos - fPositio
60
61 // Get the right child
62 // The scheme is defined by quadrant as foll
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 :
70 int bit1 = (localPosition.getY() < 0) ? 2 :
71 int bit0 = (localPosition.getZ() < 0) ? 1 :
72 return fChildren[bit0 + bit1 + bit2];
73 }
74
75 //....oooOO0OOooo........oooOO0OOooo........oo
76
77 // Non-const version of above function. Ugly b
78 OctreeNode* OctreeNode::GetChildFromPosition(c
79 {
80 G4ThreeVector localPosition = pos - fPositio
81
82 // Get the right child
83 // The scheme is defined by quadrant as foll
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 :
91 int bit1 = (localPosition.getY() < 0) ? 2 :
92 int bit0 = (localPosition.getZ() < 0) ? 1 :
93 int idx = bit0 + bit1 + bit2;
94
95 assert(idx < (int)fChildren.size());
96 return fChildren[idx];
97 }
98
99 //....oooOO0OOooo........oooOO0OOooo........oo
100
101 void OctreeNode::AddPhysicalVolume(G4VPhysical
102 {
103 // Consider adding test for the validity of
104 if (this->HasChildren()) {
105 OctreeNode* child = this->GetChildFromPosi
106 child->AddPhysicalVolume(pv);
107 }
108 else {
109 // if there are fMaxContents elements in t
110 // the bin before adding a new element.
111 if ((G4int)fContents.size() == fMaxContent
112 this->Split();
113 this->AddPhysicalVolume(pv);
114 }
115 else {
116 fContents.push_back(pv);
117 }
118 }
119 }
120
121 //....oooOO0OOooo........oooOO0OOooo........oo
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
136 fChildren[0] = new OctreeNode(newpos, hl, fM
137 newpos = G4ThreeVector(xp + xl, yp + yl, zp
138 fChildren[1] = new OctreeNode(newpos, hl, fM
139 newpos = G4ThreeVector(xp + xl, yp - yl, zp
140 fChildren[2] = new OctreeNode(newpos, hl, fM
141 newpos = G4ThreeVector(xp + xl, yp - yl, zp
142 fChildren[3] = new OctreeNode(newpos, hl, fM
143 newpos = G4ThreeVector(xp - xl, yp + yl, zp
144 fChildren[4] = new OctreeNode(newpos, hl, fM
145 newpos = G4ThreeVector(xp - xl, yp + yl, zp
146 fChildren[5] = new OctreeNode(newpos, hl, fM
147 newpos = G4ThreeVector(xp - xl, yp - yl, zp
148 fChildren[6] = new OctreeNode(newpos, hl, fM
149 newpos = G4ThreeVector(xp - xl, yp - yl, zp
150 fChildren[7] = new OctreeNode(newpos, hl, fM
151
152 // Distribute contents to children
153
154 for (int i = fContents.size() - 1; i >= 0; -
155 G4VPhysicalVolume* pv = fContents[i];
156 OctreeNode* child = this->GetChildFromPosi
157 assert(child != this);
158 child->AddPhysicalVolume(pv);
159 }
160 fContents.clear();
161 }
162
163 //....oooOO0OOooo........oooOO0OOooo........oo
164
165 std::vector<G4VPhysicalVolume*> OctreeNode::Ge
166 {
167 if (this->HasChildren()) // if has sub-node
168 {
169 std::vector<G4VPhysicalVolume*> vec;
170 std::vector<G4VPhysicalVolume*> childCont;
171 for (auto it = fChildren.begin(); it != fC
172 childCont = (*it)->GetContents();
173 for (auto jt = childCont.begin(); 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........oo
186
187 // Search for Octree nodes in a volume
188 const std::vector<G4VPhysicalVolume*> OctreeNo
189
190 {
191 // Need to search based on absolute position
192 // relative position, and maintain a list o
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()
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() - a
215 // if node within circle
216 if (d < _rad) {
217 if (aNode->HasChildren()) {
218 for (auto it = aNode->GetChildren().be
219 candidates.push_back(*it);
220 }
221 }
222 else {
223 nodes.push_back(aNode);
224 }
225 }
226 }
227
228 // G4cout << "Found " << nodes.size() << " n
229
230 // Get the physical volumes
231 G4double r2 = _rad * _rad;
232 std::vector<G4VPhysicalVolume*> pvols;
233 for (auto it = nodes.begin(); it != nodes.en
234 std::vector<G4VPhysicalVolume*> conts = (*
235 for (auto jt = conts.begin(); jt != conts.
236 if ((pos - (*jt)->GetTranslation()).mag2
237 pvols.push_back((*jt));
238 }
239 }
240 }
241 // G4cout << "Found " << pvols.size() << " c
242 return pvols;
243 }
244
245 //....oooOO0OOooo........oooOO0OOooo........oo
246
247 // Search for Octree nodes in a volume
248 void OctreeNode::SearchOctree(const G4ThreeVec
249 G4double _rad) c
250 {
251 // Need to search based on absolute position
252 // relative position, and maintain a list o
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()
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() - a
275 // if node within circle
276 if (d < _rad) {
277 if (aNode->HasChildren()) {
278 for (auto it = aNode->GetChildren().be
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.en
291 std::vector<G4VPhysicalVolume*> conts = (*
292 for (auto jt = conts.begin(); jt != conts.
293 if ((pos - (*jt)->GetTranslation()).mag2
294 out.push_back((*jt));
295 }
296 }
297 }
298 }
299
300 //....oooOO0OOooo........oooOO0OOooo........oo
301
302 const std::vector<G4VPhysicalVolume*> OctreeNo
303 {
304 const OctreeNode* child = GetChildFromPositi
305 return child->GetContents();
306 }
307
308 //....oooOO0OOooo........oooOO0OOooo........oo
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();
318 n = n + (*it)->GetNumberOfTerminalNodes();
319 }
320
321 return n;
322 }
323 //....oooOO0OOooo........oooOO0OOooo........oo
324