Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-01-18 09:58:47

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 //
0027 // Author: HoangTRAN 16/7/2019
0028 //
0029 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
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     //G4IRTUtils::GetRCutOff(1000 * CLHEP::ns) = 0.00050251
0158     //G4IRTUtils::GetRCutOff(0.1 * CLHEP::ns) = 6.4606e-06
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     //G4IRTUtils::GetRCutOff(1000 * CLHEP::ns) = 0.00050251
0243     //G4IRTUtils::GetRCutOff(0.1 * CLHEP::ns) = 6.4606e-06
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
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         //#define DEBUG
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
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 }