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