| 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 // Author: HoangTRAN 16/7/2019
28 //
29 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
30
31 template<class T,typename CONTAINER>
32 G4ThreadLocal G4OctreeFinder<T,CONTAINER> * G4OctreeFinder<T,CONTAINER>::fInstance = nullptr;
33
34 template<class T, typename CONTAINER>
35 G4OctreeFinder<T,CONTAINER> * G4OctreeFinder<T,CONTAINER>::Instance()
36 {
37 if (!fInstance)
38 {
39 fInstance = new G4OctreeFinder();
40 }
41 return fInstance;
42 }
43
44 template<class T,typename CONTAINER>
45 G4OctreeFinder<T,CONTAINER>::G4OctreeFinder()
46 : G4VFinder()
47 {}
48 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
49
50 template<class T,typename CONTAINER>
51 G4OctreeFinder<T,CONTAINER>::~G4OctreeFinder()
52 {
53 typename TreeMap::iterator it;
54 for (it = fTreeMap.begin(); it != fTreeMap.end(); it++)
55 {
56 if (it->second == nullptr)
57 {
58 continue;
59 }
60 it->second.reset();
61 }
62 fTreeMap.clear();
63 fIsOctreeBuit = false;
64 if(fInstance != nullptr)
65 {
66 delete fInstance;
67 }
68 fInstance = nullptr;
69 }
70
71 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
72
73 template<class T,typename CONTAINER>
74 void G4OctreeFinder<T,CONTAINER>::Clear()
75 {
76 typename TreeMap::iterator it;
77 for (it = fTreeMap.begin(); it != fTreeMap.end(); it++)
78 {
79 if (it->second == nullptr)
80 {
81 continue;
82 }
83 it->second.reset();
84 }
85 fTreeMap.clear();
86 fIsOctreeBuit = false;
87 }
88
89 #ifdef DEBUG
90 template<class T,typename CONTAINER>
91 void G4OctreeFinder<T,CONTAINER>::FindNaive(const G4ThreeVector& position,
92 int key,
93 G4double R,
94 std::vector<CONTAINER::iterator>& result)
95 {
96 std::map<int,PriorityList*>& listMap_test = G4ITTrackHolder::Instance()->GetLists();
97 std::map<int,PriorityList*>::iterator it = listMap_test.begin();
98 std::map<int,PriorityList*>::iterator end = listMap_test.end();
99
100 for (; it!= end; it++)
101 {
102 if(it->first == key)
103 {
104 PriorityList* Mollist=it->second;
105 result.clear();
106
107 if(Mollist == nullptr || Mollist->GetMainList() == nullptr
108 || Mollist->GetMainList()->empty() )
109 {
110 continue;
111 }
112 else
113 {
114 G4TrackList* trackList = Mollist->GetMainList();
115 G4TrackList::iterator __it = trackList->begin();
116 G4TrackList::iterator __end = trackList->end();
117 for (; __it != __end; __it++)
118 {
119 G4double r = (position - (*__it)->GetPosition()).mag();
120 if(r < R)
121 {
122 result.push_back(__it);
123 }
124 }
125 }
126 }
127 }
128 }
129 #endif
130
131 template<class T,typename CONTAINER>
132 void G4OctreeFinder<T,CONTAINER>::FindNearestInRange(const G4Track& track,
133 const G4int& key,
134 G4double R,
135 std::vector<std::pair<
136 typename CONTAINER::iterator,G4double> >&
137 result,
138 G4bool isSorted) const
139 {
140 auto it = fTreeMap.find(key);
141 if (it == fTreeMap.end())
142 {
143 return;
144 }
145 std::vector<std::pair<typename CONTAINER::iterator,G4double>> tempResult;
146 if(it->second == nullptr)
147 {
148 return;
149 }
150
151 it->second->template radiusNeighbors
152 <std::vector<std::pair<typename CONTAINER::iterator,G4double>>& >(
153 track.GetPosition(), R, tempResult);
154
155 G4int nBin = 10;
156
157 //G4IRTUtils::GetRCutOff(1000 * CLHEP::ns) = 0.00050251
158 //G4IRTUtils::GetRCutOff(0.1 * CLHEP::ns) = 6.4606e-06
159 G4double value(6.4606e-06);
160
161 G4double binOfR(std::pow(0.00050251 / value,
162 1. / static_cast<G4double>(nBin - 1)));
163
164 if( R <= 0.00050251 )
165 {
166 if(tempResult.size() < 10 && R < 0.00050251)
167 {
168 R *= binOfR;
169 #ifdef DEBUG
170 G4cout<<"recurring up R : "<< R<<" tempResult.size() : "<<tempResult.size()<<G4endl;
171 #endif
172 FindNearestInRange(track, key, R, tempResult, isSorted);
173 }
174 }
175
176 if(isSorted)
177 {
178 std::sort(tempResult.begin(),tempResult.end(),fExtractor.compareInterval);
179 }
180
181 result = std::move(tempResult);
182 }
183
184 template<class T,typename CONTAINER>
185 void G4OctreeFinder<T,CONTAINER>::FindNearest(const G4Track& track,
186 const G4int& key,
187 G4double R,
188 std::vector<std::pair<
189 typename CONTAINER::iterator,G4double> >&
190 result,
191 G4bool isSorted) const
192 {
193 auto it = fTreeMap.find(key);
194 if (it == fTreeMap.end())
195 {
196 return;
197 }
198 std::vector<std::pair<typename CONTAINER::iterator,G4double>> tempResult;
199 if(it->second == nullptr)
200 {
201 return;
202 }
203
204 it->second->template radiusNeighbors
205 <std::vector<std::pair<typename CONTAINER::iterator,G4double>>& >(
206 track.GetPosition(), R, tempResult);
207
208 if(isSorted)
209 {
210 std::sort(tempResult.begin(),tempResult.end(),fExtractor.compareInterval);
211 }
212
213 result = std::move(tempResult);
214 }
215
216
217 template<class T,typename CONTAINER>
218 void G4OctreeFinder<T,CONTAINER>::FindNearestInRange(const G4ThreeVector& position,
219 const G4int& key,
220 G4double R,
221 std::vector<std::pair<
222 typename CONTAINER::iterator,G4double> >&
223 result,
224 G4bool isSorted) const
225 {
226 typename TreeMap::const_iterator it = fTreeMap.find(key);
227 if (it == fTreeMap.end())
228 {
229 return;
230 }
231 std::vector<std::pair<typename CONTAINER::iterator,G4double>> tempResult;
232 if(it->second == nullptr)
233 {
234 return;
235 }
236
237 it->second->template radiusNeighbors
238 <std::vector<std::pair<typename CONTAINER::iterator,G4double>>& >(
239 position, R, tempResult);
240
241 G4int nBin = 10;
242 //G4IRTUtils::GetRCutOff(1000 * CLHEP::ns) = 0.00050251
243 //G4IRTUtils::GetRCutOff(0.1 * CLHEP::ns) = 6.4606e-06
244 G4double value(6.4606e-06);
245 G4double binOfR(std::pow(0.00050251 / value,
246 1. / static_cast<G4double>(nBin - 1)));
247
248 if( R <= 0.00050251 )
249 {
250 if(tempResult.size() < 10 && R < 0.00050251)
251 {
252 R *= binOfR;
253 #ifdef DEBUG
254 G4cout<<"recurring up R : "<< R<<" tempResult.size() : "<<tempResult.size()<<G4endl;
255 #endif
256 FindNearestInRange(position, key, R, tempResult, isSorted);
257 }
258 }
259
260 if(isSorted)
261 {
262 std::sort(tempResult.begin(),tempResult.end(),fExtractor.compareInterval);
263 }
264
265 result = tempResult;
266 }
267
268 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
269
270 template<class T,typename CONTAINER>
271 void G4OctreeFinder<T,CONTAINER>::
272 FindNearestInRange(const G4ThreeVector& position,
273 G4double R,
274 std::vector<std::pair<
275 typename CONTAINER::iterator,G4double> >&
276 result,
277 G4bool isSorted) const
278 {
279 typename TreeMap::const_iterator it = fTreeMap.begin();
280 std::vector<std::pair<typename CONTAINER::iterator,G4double>> Result;
281
282 for(;it != fTreeMap.end(); it++)
283 {
284 std::vector<std::pair<typename CONTAINER::iterator,G4double>> tempResult;
285 it->second->template radiusNeighbors
286 <std::vector<std::pair<typename CONTAINER::iterator,G4double>>& >(
287 position, R, tempResult);
288 if(tempResult.empty())
289 {
290 continue;
291 }
292
293 Result.insert( Result.end(), tempResult.begin(), tempResult.end() );
294
295 }
296 if(isSorted)
297 {
298 std::sort(Result.begin(),Result.end(),fExtractor.compareInterval);
299 }
300
301 result = Result;
302 }
303 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
304
305 template<class T,typename CONTAINER>
306 void G4OctreeFinder<T,CONTAINER>::BuildTreeMap(
307 const std::map<G4int,CONTAINER*>& listMap)
308 {
309 auto it = listMap.begin();
310 typename TreeMap::iterator it_Tree;
311 fTreeMap.clear();
312 for (; it!= listMap.end(); it++)
313 {
314 int key = it->first;
315 it_Tree = fTreeMap.find(key);
316 if (it_Tree != fTreeMap.end())
317 {
318 if (it_Tree->second == nullptr)
319 {
320 continue;
321 }
322 it_Tree->second.reset();
323 }
324 auto Mollist = it->second;
325 if(Mollist->empty())
326 {
327 continue;
328 }
329
330 //#define DEBUG
331 #ifdef DEBUG
332 G4cout << "** " << "Create new tree for : " << key << G4endl;
333 #endif
334 if(!Mollist->empty())
335 {
336 fTreeMap[key].reset(new Octree(Mollist->begin(),Mollist->end(),fExtractor));
337 }
338 else
339 {
340 G4ExceptionDescription exceptionDescription;
341 exceptionDescription << "should not create new tree for : " << key;
342 G4Exception("G4OCtreeFinder"
343 "::BuildReactionMap()", "BuildReactionMap02",
344 FatalException, exceptionDescription);
345 }
346
347 #ifdef DEBUG
348
349 auto __it = Mollist->begin();
350 auto __end = Mollist->end();
351
352 for (; __it != __end; __it++)
353 {
354 G4cout<< "molecule : "<<(*__it)->GetPosition()<< G4endl;
355 }
356 #endif
357 #undef DEBUG
358
359 }
360 fIsOctreeBuit = true;
361 }
362
363 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
364
365 template<class T,typename CONTAINER>
366 void G4OctreeFinder<T,CONTAINER>::SetOctreeUsed(G4bool used)
367 {
368 fIsOctreeUsed = used;
369 }
370
371 template<class T,typename CONTAINER>
372 G4bool G4OctreeFinder<T,CONTAINER>::IsOctreeUsed() const
373 {
374 return fIsOctreeUsed;
375 }
376
377 template<class T,typename CONTAINER>
378 void G4OctreeFinder<T,CONTAINER>::SetOctreeBuilt(G4bool used)
379 {
380 fIsOctreeBuit = used;
381 }
382
383 template<class T,typename CONTAINER>
384 G4bool G4OctreeFinder<T,CONTAINER>::IsOctreeBuilt() const
385 {
386 return fIsOctreeBuit;
387 }
388