File indexing completed on 2025-01-18 09:58:47
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031 template<class T,typename CONTAINER>
0032 G4ThreadLocal G4OctreeFinder<T,CONTAINER> * G4OctreeFinder<T,CONTAINER>::fInstance = nullptr;
0033
0034 template<class T, typename CONTAINER>
0035 G4OctreeFinder<T,CONTAINER> * G4OctreeFinder<T,CONTAINER>::Instance()
0036 {
0037 if (!fInstance)
0038 {
0039 fInstance = new G4OctreeFinder();
0040 }
0041 return fInstance;
0042 }
0043
0044 template<class T,typename CONTAINER>
0045 G4OctreeFinder<T,CONTAINER>::G4OctreeFinder()
0046 : G4VFinder()
0047 {}
0048
0049
0050 template<class T,typename CONTAINER>
0051 G4OctreeFinder<T,CONTAINER>::~G4OctreeFinder()
0052 {
0053 typename TreeMap::iterator it;
0054 for (it = fTreeMap.begin(); it != fTreeMap.end(); it++)
0055 {
0056 if (it->second == nullptr)
0057 {
0058 continue;
0059 }
0060 it->second.reset();
0061 }
0062 fTreeMap.clear();
0063 fIsOctreeBuit = false;
0064 if(fInstance != nullptr)
0065 {
0066 delete fInstance;
0067 }
0068 fInstance = nullptr;
0069 }
0070
0071
0072
0073 template<class T,typename CONTAINER>
0074 void G4OctreeFinder<T,CONTAINER>::Clear()
0075 {
0076 typename TreeMap::iterator it;
0077 for (it = fTreeMap.begin(); it != fTreeMap.end(); it++)
0078 {
0079 if (it->second == nullptr)
0080 {
0081 continue;
0082 }
0083 it->second.reset();
0084 }
0085 fTreeMap.clear();
0086 fIsOctreeBuit = false;
0087 }
0088
0089 #ifdef DEBUG
0090 template<class T,typename CONTAINER>
0091 void G4OctreeFinder<T,CONTAINER>::FindNaive(const G4ThreeVector& position,
0092 int key,
0093 G4double R,
0094 std::vector<CONTAINER::iterator>& result)
0095 {
0096 std::map<int,PriorityList*>& listMap_test = G4ITTrackHolder::Instance()->GetLists();
0097 std::map<int,PriorityList*>::iterator it = listMap_test.begin();
0098 std::map<int,PriorityList*>::iterator end = listMap_test.end();
0099
0100 for (; it!= end; it++)
0101 {
0102 if(it->first == key)
0103 {
0104 PriorityList* Mollist=it->second;
0105 result.clear();
0106
0107 if(Mollist == nullptr || Mollist->GetMainList() == nullptr
0108 || Mollist->GetMainList()->empty() )
0109 {
0110 continue;
0111 }
0112 else
0113 {
0114 G4TrackList* trackList = Mollist->GetMainList();
0115 G4TrackList::iterator __it = trackList->begin();
0116 G4TrackList::iterator __end = trackList->end();
0117 for (; __it != __end; __it++)
0118 {
0119 G4double r = (position - (*__it)->GetPosition()).mag();
0120 if(r < R)
0121 {
0122 result.push_back(__it);
0123 }
0124 }
0125 }
0126 }
0127 }
0128 }
0129 #endif
0130
0131 template<class T,typename CONTAINER>
0132 void G4OctreeFinder<T,CONTAINER>::FindNearestInRange(const G4Track& track,
0133 const G4int& key,
0134 G4double R,
0135 std::vector<std::pair<
0136 typename CONTAINER::iterator,G4double> >&
0137 result,
0138 G4bool isSorted) const
0139 {
0140 auto it = fTreeMap.find(key);
0141 if (it == fTreeMap.end())
0142 {
0143 return;
0144 }
0145 std::vector<std::pair<typename CONTAINER::iterator,G4double>> tempResult;
0146 if(it->second == nullptr)
0147 {
0148 return;
0149 }
0150
0151 it->second->template radiusNeighbors
0152 <std::vector<std::pair<typename CONTAINER::iterator,G4double>>& >(
0153 track.GetPosition(), R, tempResult);
0154
0155 G4int nBin = 10;
0156
0157
0158
0159 G4double value(6.4606e-06);
0160
0161 G4double binOfR(std::pow(0.00050251 / value,
0162 1. / static_cast<G4double>(nBin - 1)));
0163
0164 if( R <= 0.00050251 )
0165 {
0166 if(tempResult.size() < 10 && R < 0.00050251)
0167 {
0168 R *= binOfR;
0169 #ifdef DEBUG
0170 G4cout<<"recurring up R : "<< R<<" tempResult.size() : "<<tempResult.size()<<G4endl;
0171 #endif
0172 FindNearestInRange(track, key, R, tempResult, isSorted);
0173 }
0174 }
0175
0176 if(isSorted)
0177 {
0178 std::sort(tempResult.begin(),tempResult.end(),fExtractor.compareInterval);
0179 }
0180
0181 result = tempResult;
0182 }
0183
0184 template<class T,typename CONTAINER>
0185 void G4OctreeFinder<T,CONTAINER>::FindNearest(const G4Track& track,
0186 const G4int& key,
0187 G4double R,
0188 std::vector<std::pair<
0189 typename CONTAINER::iterator,G4double> >&
0190 result,
0191 G4bool isSorted) const
0192 {
0193 auto it = fTreeMap.find(key);
0194 if (it == fTreeMap.end())
0195 {
0196 return;
0197 }
0198 std::vector<std::pair<typename CONTAINER::iterator,G4double>> tempResult;
0199 if(it->second == nullptr)
0200 {
0201 return;
0202 }
0203
0204 it->second->template radiusNeighbors
0205 <std::vector<std::pair<typename CONTAINER::iterator,G4double>>& >(
0206 track.GetPosition(), R, tempResult);
0207
0208 if(isSorted)
0209 {
0210 std::sort(tempResult.begin(),tempResult.end(),fExtractor.compareInterval);
0211 }
0212
0213 result = tempResult;
0214 }
0215
0216
0217 template<class T,typename CONTAINER>
0218 void G4OctreeFinder<T,CONTAINER>::FindNearestInRange(const G4ThreeVector& position,
0219 const G4int& key,
0220 G4double R,
0221 std::vector<std::pair<
0222 typename CONTAINER::iterator,G4double> >&
0223 result,
0224 G4bool isSorted) const
0225 {
0226 typename TreeMap::const_iterator it = fTreeMap.find(key);
0227 if (it == fTreeMap.end())
0228 {
0229 return;
0230 }
0231 std::vector<std::pair<typename CONTAINER::iterator,G4double>> tempResult;
0232 if(it->second == nullptr)
0233 {
0234 return;
0235 }
0236
0237 it->second->template radiusNeighbors
0238 <std::vector<std::pair<typename CONTAINER::iterator,G4double>>& >(
0239 position, R, tempResult);
0240
0241 G4int nBin = 10;
0242
0243
0244 G4double value(6.4606e-06);
0245 G4double binOfR(std::pow(0.00050251 / value,
0246 1. / static_cast<G4double>(nBin - 1)));
0247
0248 if( R <= 0.00050251 )
0249 {
0250 if(tempResult.size() < 10 && R < 0.00050251)
0251 {
0252 R *= binOfR;
0253 #ifdef DEBUG
0254 G4cout<<"recurring up R : "<< R<<" tempResult.size() : "<<tempResult.size()<<G4endl;
0255 #endif
0256 FindNearestInRange(position, key, R, tempResult, isSorted);
0257 }
0258 }
0259
0260 if(isSorted)
0261 {
0262 std::sort(tempResult.begin(),tempResult.end(),fExtractor.compareInterval);
0263 }
0264
0265 result = tempResult;
0266 }
0267
0268
0269
0270 template<class T,typename CONTAINER>
0271 void G4OctreeFinder<T,CONTAINER>::
0272 FindNearestInRange(const G4ThreeVector& position,
0273 G4double R,
0274 std::vector<std::pair<
0275 typename CONTAINER::iterator,G4double> >&
0276 result,
0277 G4bool isSorted) const
0278 {
0279 typename TreeMap::const_iterator it = fTreeMap.begin();
0280 std::vector<std::pair<typename CONTAINER::iterator,G4double>> Result;
0281
0282 for(;it != fTreeMap.end(); it++)
0283 {
0284 std::vector<std::pair<typename CONTAINER::iterator,G4double>> tempResult;
0285 it->second->template radiusNeighbors
0286 <std::vector<std::pair<typename CONTAINER::iterator,G4double>>& >(
0287 position, R, tempResult);
0288 if(tempResult.empty())
0289 {
0290 continue;
0291 }
0292
0293 Result.insert( Result.end(), tempResult.begin(), tempResult.end() );
0294
0295 }
0296 if(isSorted)
0297 {
0298 std::sort(Result.begin(),Result.end(),fExtractor.compareInterval);
0299 }
0300
0301 result = Result;
0302 }
0303
0304
0305 template<class T,typename CONTAINER>
0306 void G4OctreeFinder<T,CONTAINER>::BuildTreeMap(
0307 const std::map<G4int,CONTAINER*>& listMap)
0308 {
0309 auto it = listMap.begin();
0310 typename TreeMap::iterator it_Tree;
0311 fTreeMap.clear();
0312 for (; it!= listMap.end(); it++)
0313 {
0314 int key = it->first;
0315 it_Tree = fTreeMap.find(key);
0316 if (it_Tree != fTreeMap.end())
0317 {
0318 if (it_Tree->second == nullptr)
0319 {
0320 continue;
0321 }
0322 it_Tree->second.reset();
0323 }
0324 auto Mollist = it->second;
0325 if(Mollist->empty())
0326 {
0327 continue;
0328 }
0329
0330
0331 #ifdef DEBUG
0332 G4cout << "** " << "Create new tree for : " << key << G4endl;
0333 #endif
0334 if(!Mollist->empty())
0335 {
0336 fTreeMap[key].reset(new Octree(Mollist->begin(),Mollist->end(),fExtractor));
0337 }
0338 else
0339 {
0340 G4ExceptionDescription exceptionDescription;
0341 exceptionDescription << "should not create new tree for : " << key;
0342 G4Exception("G4OCtreeFinder"
0343 "::BuildReactionMap()", "BuildReactionMap02",
0344 FatalException, exceptionDescription);
0345 }
0346
0347 #ifdef DEBUG
0348
0349 auto __it = Mollist->begin();
0350 auto __end = Mollist->end();
0351
0352 for (; __it != __end; __it++)
0353 {
0354 G4cout<< "molecule : "<<(*__it)->GetPosition()<< G4endl;
0355 }
0356 #endif
0357 #undef DEBUG
0358
0359 }
0360 fIsOctreeBuit = true;
0361 }
0362
0363
0364
0365 template<class T,typename CONTAINER>
0366 void G4OctreeFinder<T,CONTAINER>::SetOctreeUsed(G4bool used)
0367 {
0368 fIsOctreeUsed = used;
0369 }
0370
0371 template<class T,typename CONTAINER>
0372 G4bool G4OctreeFinder<T,CONTAINER>::IsOctreeUsed() const
0373 {
0374 return fIsOctreeUsed;
0375 }
0376
0377 template<class T,typename CONTAINER>
0378 void G4OctreeFinder<T,CONTAINER>::SetOctreeBuilt(G4bool used)
0379 {
0380 fIsOctreeBuit = used;
0381 }
0382
0383 template<class T,typename CONTAINER>
0384 G4bool G4OctreeFinder<T,CONTAINER>::IsOctreeBuilt() const
0385 {
0386 return fIsOctreeBuit;
0387 }