CMakeLists.txt

@@ -464,7 +464,8 @@ if (BENCHMARK)
   set(SRC_CPP ${SRC_CPP}
     source/benchmarking/BenchmarkResult.cpp
     source/benchmarking/Benchmarker.cpp
-    #source/benchmarking/NavigationBenchmarker.cpp
+    source/benchmarking/NavigationBenchmarker.cpp
+    source/benchmarking/NavigationBenchmarker.cpp
     source/benchmarking/VolumePointers.cpp
   )
   set(SRC_CUDA ${SRC_CUDA}
@@ -564,7 +565,7 @@ set(TEST_EXECUTABLES_ROOT
 
   # higher level benchmarks
   #${CMAKE_SOURCE_DIR}/test/globalbenchmarks/LocatePointsBenchmark.cpp
-  #${CMAKE_SOURCE_DIR}/test/globalbenchmarks/NavigationBenchmark.cpp
+  ${CMAKE_SOURCE_DIR}/test/globalbenchmarks/NavigationBenchmark.cpp
 )
 if(BENCHMARK)
   set(TEST_EXECUTABLES_ROOT

VecCore/inc/BitSet.h

 /// \file VecCore/BitSet.h
 /// \author Philippe Canal (pcanal@fnal.gov)
-/// \author Exported from the original verison in ROOT (root.cern.ch).
+/// \author Exported from the original version in ROOT (root.cern.ch).
 
 #ifndef VECCORE_BITSET_H
 #define VECCORE_BITSET_H
@@ -73,7 +73,7 @@ namespace VecCore {
       BitSet(size_t new_size, const BitSet &other) : fNbits(other.fNbits), fData(new_size, other.fData) {
          // We assume that the memory allocated and use is large enough to
          // hold the full values (i.e. Sizeof(new_size) )
-         // If new_size is smaller than the size of 'other' then the copy is trucated.
+         // If new_size is smaller than the size of 'other' then the copy is truncated.
          // if new_size is greater than the size of 'other' the value are zero initialized.
 
          if (fNbits*8 > new_size) {

VecCore/inc/VariableSizeObj.h

@@ -86,9 +86,7 @@ namespace VecCore {
       VECGEOM_CUDA_HEADER_BOTH
       static  Cont *MakeInstance(size_t nvalues, const T&... params) {
          // Make an instance of the class which allocates the node array. To be
-         // released using ReleaseInstance. If addr is non-zero, the user promised that
-         // addr contains at least that many bytes:  size_t needed = SizeOf(nvalues);
-
+         // released using ReleaseInstance.
          size_t needed = SizeOf(nvalues);
          char *ptr = new char[ needed ];
          if (!ptr) return 0;
@@ -114,7 +112,7 @@ namespace VecCore {
       // The equivalent of the copy constructor
       static Cont *MakeCopy(const Cont &other)
       {
-         // Make a copy of a the variable size array and its container.
+         // Make a copy of the variable size array and its container.
 
          size_t needed = SizeOf(other.GetVariableData().fN);
          char *ptr = new char[ needed ];

benchmarking/NavigationBenchmarker.h

@@ -45,9 +45,7 @@ namespace vecgeom {
     SOA3D<Precision> const& dirs
     );
 
-  void testVectorSafety( VPlacedVolume const* top );
-
-  void runNavigationBenchmarks( VPlacedVolume const* top, int np, int nreps, Precision bias = 0.8);
+  void runNavigationBenchmarks( LogicalVolume const* top, int np, int nreps, Precision bias = 0.8);
 
   bool validateNavigationStepAgainstRoot(
     Vector3D<Precision> const& pos,

navigation/NavigationState.h

@@ -106,21 +106,21 @@ private:
 
 
 public:
-    // replaces the volume pointers from CPU volumes in fPath
-     // to the equivalent pointers on the GPU
-     // uses the CudaManager to do so
-    void ConvertToGPUPointers();
-
-    // replaces the pointers from GPU volumes in fPath
-    // to the equivalent pointers on the CPU
-    // uses the CudaManager to do so
-    void ConvertToCPUPointers();
-
-  // Enumerate the part of the private interface, we want to expose.
-  using Base_t::MakeCopy;
-  using Base_t::MakeCopyAt;
-  using Base_t::ReleaseInstance;
-  using Base_t::SizeOf;
+   // replaces the volume pointers from CPU volumes in fPath
+   // to the equivalent pointers on the GPU
+   // uses the CudaManager to do so
+   void ConvertToGPUPointers();
+
+   // replaces the pointers from GPU volumes in fPath
+   // to the equivalent pointers on the CPU
+   // uses the CudaManager to do so
+   void ConvertToCPUPointers();
+
+   // Enumerate the part of the private interface, we want to expose.
+   using Base_t::MakeCopy;
+   using Base_t::MakeCopyAt;
+   using Base_t::ReleaseInstance;
+   using Base_t::SizeOf;
 
    // produces a compact navigation state object of a certain depth
    // the caller can give a memory address where the object will
@@ -128,7 +128,7 @@ public:
    // the caller has to make sure that the size of the external memory
    // is >= sizeof(NavigationState) + sizeof(VPlacedVolume*)*maxlevel
    //
-   // Both MakeInstance, MakeInstanceAt, MakeCopy and MakeCopyAT are provided by 
+   // Methods MakeInstance(), MakeInstanceAt(), MakeCopy() and MakeCopyAt() are provided by
    // VariableSizeObjectInterface
 
    VECGEOM_CUDA_HEADER_BOTH
@@ -432,32 +432,17 @@ VECGEOM_INLINE
 VECGEOM_CUDA_HEADER_BOTH
 void NavigationState::Print() const
 {
-   // printf("bool: fOnBoundary=%i %p (%l bytes)\n", fOnBoundary, static_cast<void*>(fOnBoundary), sizeof(bool));
-   // printf("Transf3D: matrix (%l bytes)\n", sizeof(Transformation3D) );
    // printf("VariableSizeObj: fPath=%p (%l bytes)\n", fPath, sizeof(fPath));
-
 #ifndef VECGEOM_NVCC
    printf("NavState: Level(cur/max)=%i/%i,  onBoundary=%s, topVol=<%s>, this=%p\n",
-          fCurrentLevel,
-          GetMaxLevel(),
-          (fOnBoundary?"true":"false"),
-          (Top()?
-           Top()->GetLabel().c_str():
-           "NULL"),
-          this );
+          fCurrentLevel, GetMaxLevel(), (fOnBoundary?"true":"false"),
+          (Top()? Top()->GetLabel().c_str():"NULL"), this );
 #else
    printf("NavState: Level(cur/max)=%i/%i,  onBoundary=%s, topVol=<%p>, this=%p\n",
-      fCurrentLevel,
-      GetMaxLevel(),
-      (fOnBoundary?"true":"false"),
-      Top(),
-      this );
+          fCurrentLevel, GetMaxLevel(), (fOnBoundary?"true":"false"),
+          Top(), this );
 #endif
 
-   // std::cerr << "NavState: Level(cur/max)=" << fCurrentLevel <<'/'<< GetMaxLevel()
-   //           <<" onBoundary="<< fOnBoundary
-   //           <<" topVol="<< Top() <<" this="<< this
-   //           << std::endl;
 }
 
 

navigation/SimpleNavigator.h

@@ -12,6 +12,7 @@
 #include "base/Vector3D.h"
 #include "management/GeoManager.h"
 #include "navigation/NavigationState.h"
+#include "navigation/NavStatePool.h"
 
 #ifdef VECGEOM_ROOT
 #include "management/RootGeoManager.h"
@@ -136,8 +137,8 @@ public:
          Container3D const & /*global dirs*/,
          Container3D & /*workspace for localpoints*/,
          Container3D & /*workspace for localdirs*/,
-         NavigationState **  /* array of pointers to NavigationStates for currentstates */,
-         NavigationState **  /* array of pointers to NabigationStates for outputstates */,
+         NavStatePool const&  /* array of pointers to NavigationStates for currentstates */,
+         NavStatePool &  /* array of pointers to NabigationStates for outputstates */,
          Precision const * /* pSteps -- proposed steps */,
          Precision * /* safeties */,
          Precision * /* distances; steps */,
@@ -160,8 +161,8 @@ public:
     * and a global point; mainly for debugging purposes
     */
    void InspectSafetyForPoint(
-		   Vector3D<Precision> const & /* global point */,
-		   NavigationState const & /* current state */
+           Vector3D<Precision> const & /* global point */,
+           NavigationState const & /* current state */
       ) const;
 
 }; // end of class declaration
@@ -400,15 +401,14 @@ void SimpleNavigator::GetSafeties(Container3D const & globalpoints,
  * Navigation interface for baskets; templates on Container3D which might be a SOA3D or AOS3D container
  */
 
-
 template <typename Container3D>
 void SimpleNavigator::FindNextBoundaryAndStep(
          Container3D const & globalpoints,
          Container3D const & globaldirs,
          Container3D       & localpoints,
          Container3D       & localdirs,
-         NavigationState  ** currentstates,
-         NavigationState  ** newstates,
+         NavStatePool const& currentstates,
+         NavStatePool      & newstates,
          Precision const   * pSteps,
          Precision         * safeties,
          Precision         * distances,

source/CudaManager.cpp

@@ -194,7 +194,7 @@ void CudaManager::CleanGpu() {
 template <typename Coll>
 bool CudaManager::AllocateCollectionOnCoproc(const char *verbose_title,
                                              const Coll &data,
-					     bool isforplacedvol
+                         bool isforplacedvol
                                              )
 {
    // NOTE: Code need to be enhanced to propage the error correctly.
@@ -213,7 +213,7 @@ bool CudaManager::AllocateCollectionOnCoproc(const char *verbose_title,
    for (auto i : data) {
       memory_map[ToCpuAddress(i)] = gpu_address;
       if(isforplacedvol)
-	  fGPUtoCPUmapForPlacedVolumes[ gpu_address ] = i;
+      fGPUtoCPUmapForPlacedVolumes[ gpu_address ] = i;
       gpu_address += i->DeviceSizeOf();
    }
 

source/GeoManager.cpp

@@ -48,6 +48,7 @@ VPlacedVolume* GeoManager::FindPlacedVolume(const int id) {
 VPlacedVolume* GeoManager::FindPlacedVolume(char const *const label) {
   VPlacedVolume *output = NULL;
   bool multiple = false;
+  printf("FindPlacedVolume: %s - fPlacedVolumesMap.size() = %li\n", label, fPlacedVolumesMap.size());
   for (auto v = fPlacedVolumesMap.begin(), v_end = fPlacedVolumesMap.end();
        v != v_end; ++v) {
     if (v->second->GetLabel() == label) {
@@ -79,7 +80,12 @@ LogicalVolume* GeoManager::FindLogicalVolume(char const *const label) {
   bool multiple = false;
   for (auto v = fLogicalVolumesMap.begin(), v_end = fLogicalVolumesMap.end();
        v != v_end; ++v) {
-    if (v->second->GetLabel() == label) {
+
+    // This is needed for CMS... strip off last 4 chars from logVol name
+    const std::string& fullname = v->second->GetLabel();
+    std::string strippedName( fullname, 0, fullname.length()-4 );
+
+    if (strippedName.compare(label)==0) {
       if (!output) {
         output = v->second;
       } else {

source/benchmarking/Benchmarker.cpp

@@ -6,7 +6,6 @@
 #include "base/SOA3D.h"
 #include "base/Stopwatch.h"
 #include "base/Transformation3D.h"
-#include "volumes/LogicalVolume.h"
 #include "volumes/PlacedBox.h"
 #include "volumes/utilities/VolumeUtilities.h"
 

source/benchmarking/NavigationBenchmarker.cpp

@@ -44,8 +44,7 @@ namespace vecgeom {
 
 //==================================
 
-Precision benchmarkLocatePoint( VPlacedVolume_t startVol, int nPoints, int nReps,
-                                SOA3D<Precision> const& points) {
+Precision benchmarkLocatePoint( int nPoints, int nReps, SOA3D<Precision> const& points) {
 
   NavigationState * state = NavigationState::MakeInstance( GeoManager::Instance().getMaxDepth() );
   vecgeom::SimpleNavigator nav;
@@ -55,7 +54,7 @@ Precision benchmarkLocatePoint( VPlacedVolume_t startVol, int nPoints, int nReps
   for(int n=0; n<nReps; ++n) {
     for( int i=0; i<nPoints; ++i ) {
       state->Clear();
-      nav.LocatePoint( startVol, points[i], *state, startVol==GeoManager::Instance().GetWorld());
+      nav.LocatePoint( GeoManager::Instance().GetWorld(), points[i], *state, true);
     }
   }
   Precision elapsed = timer.Stop();
@@ -66,7 +65,7 @@ Precision benchmarkLocatePoint( VPlacedVolume_t startVol, int nPoints, int nReps
 
 //==================================
 
-Precision benchmarkSerialNavigation( VPlacedVolume_t startVol, int nPoints, int nReps,
+Precision benchmarkSerialNavigation( int nPoints, int nReps,
                                      SOA3D<Precision> const& points,
                                      SOA3D<Precision> const& dirs ) {
 
@@ -80,7 +79,7 @@ Precision benchmarkSerialNavigation( VPlacedVolume_t startVol, int nPoints, int
 
   for( int i=0; i<nPoints; ++i ) {
     curStates[i]->Clear();
-    nav.LocatePoint( startVol, points[i], *curStates[i], startVol==GeoManager::Instance().GetWorld());
+    nav.LocatePoint( GeoManager::Instance().GetWorld(), points[i], *curStates[i], true);
   }
 
   Precision maxStep = kInfinity, step=0.0;
@@ -107,14 +106,10 @@ Precision benchmarkSerialNavigation( VPlacedVolume_t startVol, int nPoints, int
 
 //==================================
 
-Precision benchmarkVectorNavigation( VPlacedVolume_t startVol, int nPoints, int nReps,
+Precision benchmarkVectorNavigation( int nPoints, int nReps,
                                      SOA3D<Precision> const& points,
                                      SOA3D<Precision> const& dirs ) {
 
-  // now setup all the navigation states
-  NavigationState ** curStates = new NavigationState*[nPoints];
-  NavigationState ** newStates = new NavigationState*[nPoints];
-
   SOA3D<Precision> workspace1(nPoints);
   SOA3D<Precision> workspace2(nPoints);
 
@@ -128,10 +123,17 @@ Precision benchmarkVectorNavigation( VPlacedVolume_t startVol, int nPoints, int
   memset(safeties, 0, sizeof(Precision)*nPoints);
 
   vecgeom::SimpleNavigator nav;
-  for (int i=0;i<nPoints;++i){
-    curStates[i] = NavigationState::MakeInstance( GeoManager::Instance().getMaxDepth() );
-    newStates[i] = NavigationState::MakeInstance( GeoManager::Instance().getMaxDepth() );
-    nav.LocatePoint( startVol, points[i], *curStates[i], startVol==GeoManager::Instance().GetWorld());
+
+  // setup all the navigation states
+  NavStatePool curStates(nPoints, GeoManager::Instance().getMaxDepth() );
+  NavStatePool newStates(nPoints, GeoManager::Instance().getMaxDepth() );
+  //NavigationState ** curStates = new NavigationState*[nPoints];
+  //NavigationState ** newStates = new NavigationState*[nPoints];
+  for( int i=0; i<nPoints; ++i ) {
+    curStates[i]->Clear();
+    // curStates[i] = NavigationState::MakeInstance( GeoManager::Instance().getMaxDepth() );
+    // newStates[i] = NavigationState::MakeInstance( GeoManager::Instance().getMaxDepth() );
+    nav.LocatePoint( GeoManager::Instance().GetWorld(), points[i], *curStates[i], true);
   }
 
   Stopwatch timer;
@@ -147,8 +149,8 @@ Precision benchmarkVectorNavigation( VPlacedVolume_t startVol, int nPoints, int
     NavigationState::ReleaseInstance( curStates[i] );
     NavigationState::ReleaseInstance( newStates[i] );
   }
-  delete[] curStates;
-  delete[] newStates;
+  // delete[] curStates;
+  // delete[] newStates;
 
   _mm_free(intworkspace);
   _mm_free(maxSteps);
@@ -240,58 +242,25 @@ Precision benchmarkGeant4Navigation( int nPoints, int nReps,
 }
 #endif
 
-//==================================
-// function to test safety
-void testVectorSafety( VPlacedVolume_t startVol ){
-   SOA3D<Precision> points(1024);
-   SOA3D<Precision> workspace(1024);
-   Precision * safeties = (Precision *) _mm_malloc(sizeof(Precision)*1024,32);
-   vecgeom::volumeUtilities::FillUncontainedPoints( *startVol, points );
-
-   // now setup all the navigation states
-   NavigationState ** states = new NavigationState*[1024];
-
-   vecgeom::SimpleNavigator nav;
-
-   for (int i=0;i<1024;++i){
-     states[i] = NavigationState::MakeInstance( GeoManager::Instance().getMaxDepth() );
-     nav.LocatePoint( startVol, points[i], *states[i], startVol==GeoManager::Instance().GetWorld());
-   }
-
-   // calculate safeties with vector interface
-   nav.GetSafeties(points, states, workspace, safeties );
-
-   // verify against serial interface
-   for (int i=0;i<1024;++i){
-       vecgeom::Assert( safeties[i] == nav.GetSafety( points[i], *states[i] ), ""
-                        " Problem in VectorSafety (in SimpleNavigator)" );
-   }
-   std::cout << "Safety test passed\n";
-
-   for (int i=0;i<1024;++i) NavigationState::ReleaseInstance( states[i] );
-   delete states;
-   _mm_free(safeties);
-}
-
 //=======================================
 /// Function to run navigation benchmarks
-void runNavigationBenchmarks( VPlacedVolume_t startVol, int np, int nreps, Precision bias) {
+void runNavigationBenchmarks( LogicalVolume const* startVol, int np, int nreps, Precision bias) {
 
   SOA3D<Precision> points(np);
   SOA3D<Precision> locpts(np);
   SOA3D<Precision> dirs(np);
-  vecgeom::volumeUtilities::FillGlobalPointsAndDirectionsForLogicalVolume(
-    startVol->GetLogicalVolume(), locpts, points, dirs, bias, np);
+
+  vecgeom::volumeUtilities::FillGlobalPointsAndDirectionsForLogicalVolume( startVol, locpts, points, dirs, bias, np);
 
   Precision cputime;
 
-  cputime = benchmarkLocatePoint(startVol,np,nreps,points);
+  cputime = benchmarkLocatePoint(np,nreps,points);
   printf("CPU elapsed time (locating and setting steps) %f ms\n", 1000.*cputime);
 
-  cputime = benchmarkSerialNavigation(startVol,np,nreps,points, dirs);
+  cputime = benchmarkSerialNavigation(np,nreps,points, dirs);
   printf("CPU elapsed time (serialized navigation) %f ms\n", 1000.*cputime);
 
-  cputime = benchmarkVectorNavigation(startVol,np,nreps,points, dirs);
+  cputime = benchmarkVectorNavigation(np,nreps,points, dirs);
   printf("CPU elapsed time (vectorized navigation) %f ms\n", 1000.*cputime);
 
 #ifdef VECGEOM_ROOT
@@ -336,6 +305,7 @@ bool validateNavigationStepAgainstRoot(
       std::cerr << " OUTSIDEERROR \n";
     }
   }
+
   else if( Abs(testStep - rootnav->GetStep()) > 100.*kTolerance
            || rootnav->GetCurrentNode() != RootGeoManager::Instance().tgeonode(testState.Top()) )
   {
@@ -428,11 +398,12 @@ bool validateVecGeomNavigation( int np, SOA3D<Precision> const& points, SOA3D<Pr
   bool result = true;
 
   // now setup all the navigation states - one loop at a time for better data locality
-  NavigationState** origStates     = new NavigationState*[np];
-  for( int i=0; i<np; ++i) origStates[i] = NavigationState::MakeInstance( GeoManager::Instance().getMaxDepth() );
-
-  NavigationState** vgSerialStates = new NavigationState*[np];
-  for( int i=0; i<np; ++i) vgSerialStates[i] = NavigationState::MakeInstance( GeoManager::Instance().getMaxDepth() );
+  NavStatePool origStates(np, GeoManager::Instance().getMaxDepth() );
+  NavStatePool vgSerialStates(np, GeoManager::Instance().getMaxDepth() );
+  // NavigationState** origStates     = new NavigationState*[np];
+  // for( int i=0; i<np; ++i) origStates[i] = NavigationState::MakeInstance( GeoManager::Instance().getMaxDepth() );
+  // NavigationState** vgSerialStates = new NavigationState*[np];
+  // for( int i=0; i<np; ++i) vgSerialStates[i] = NavigationState::MakeInstance( GeoManager::Instance().getMaxDepth() );
 
   vecgeom::SimpleNavigator nav;
   Precision * maxSteps = (Precision *) _mm_malloc(sizeof(Precision)*np,32);
@@ -473,9 +444,6 @@ bool validateVecGeomNavigation( int np, SOA3D<Precision> const& points, SOA3D<Pr
 
     result &= ok;
     if( !ok ) {
-
-      TGeoNavigator* rootnav = ::gGeoManager->GetCurrentNavigator();
-
       std::cout << "\n=======> Summary: ITERATION " << i
                 <<" - pos = " << pos <<" dir = " << dir
                 <<" / Steps (";
@@ -488,6 +456,7 @@ bool validateVecGeomNavigation( int np, SOA3D<Precision> const& points, SOA3D<Pr
       std::cout<<"VecGeom): ";
 
 #ifdef VECGEOM_ROOT
+      TGeoNavigator* rootnav = ::gGeoManager->GetCurrentNavigator();
       std::cout<< rootnav->GetStep() << " / ";
 #endif
 #ifdef VECGEOM_GEANT4
@@ -504,6 +473,7 @@ bool validateVecGeomNavigation( int np, SOA3D<Precision> const& points, SOA3D<Pr
       std::cout << (nextPV ? nextPV->GetName() : "NULL") <<" / ";
 #endif
       std::cout << (vgSerialStates[i]->Top()? vgSerialStates[i]->Top()->GetLabel() : "NULL") << "\n";
+
       // nav.InspectEnvironmentForPointAndDirection( pos, dir, *origState );
     }
   }
@@ -517,8 +487,9 @@ bool validateVecGeomNavigation( int np, SOA3D<Precision> const& points, SOA3D<Pr
   //=== N-particle navigation interface
 
   //--- Creating vgVectorStates
-  NavigationState** vgVectorStates = new NavigationState*[np];
-  for( int i=0; i<np; ++i) vgVectorStates[i] = NavigationState::MakeInstance( GeoManager::Instance().getMaxDepth() );
+  NavStatePool vgVectorStates(np, GeoManager::Instance().getMaxDepth() );
+  // NavigationState** vgVectorStates = new NavigationState*[np];
+  // for( int i=0; i<np; ++i) vgVectorStates[i] = NavigationState::MakeInstance( GeoManager::Instance().getMaxDepth() );
 
   SOA3D<Precision> workspace1(np);
   SOA3D<Precision> workspace2(np);

test/benchmark/BoxBenchmark.cpp

@@ -8,7 +8,7 @@ using namespace vecgeom;
 
 int main(int argc, char* argv[]) {
   OPTION_INT(npoints,1024);
-  OPTION_INT(nrep,1024);
+  OPTION_INT(nrep,4);
   OPTION_DOUBLE(dx,1.);
   OPTION_DOUBLE(dy,2.);
   OPTION_DOUBLE(dz,3.);

test/benchmark/ParaboloidBenchmark.cpp

@@ -25,8 +25,8 @@ int main() {
     
     Benchmarker tester(GeoManager::Instance().GetWorld());
     tester.SetVerbosity(3);
-    tester.SetPointCount(1<<13);
-    tester.SetRepetitions(1024);
+    tester.SetPointCount(1<<10);
+    tester.SetRepetitions(4);
     std::cout<<"Prepared to run benchmarker\n";
     tester.RunBenchmark();
     

test/benchmark/PolyconeBenchmark.cpp

@@ -16,7 +16,7 @@ using namespace vecgeom;
 
 int main(int argc, char* argv[]) {
   OPTION_INT(npoints,1024);
-  OPTION_INT(nrep,1024);
+  OPTION_INT(nrep,4);
   OPTION_DOUBLE(phistart,0.);
   OPTION_DOUBLE(phidelta,kTwoPi);
 

test/benchmark/PolyhedronBenchmark.cpp

@@ -58,8 +58,8 @@ int main(int nArgs, char **args) {
     Benchmarker benchmarker(world);
     benchmarker.SetVerbosity(3);
     benchmarker.SetPoolMultiplier(1);
-    benchmarker.SetRepetitions(8192);
-    benchmarker.SetPointCount(128);
+    benchmarker.SetRepetitions(4);
+    benchmarker.SetPointCount(1024);
     benchmarker.RunInsideBenchmark();
     benchmarker.RunToInBenchmark();
     benchmarker.RunToOutBenchmark();

test/benchmark/TorusBenchmark.cpp

@@ -10,7 +10,7 @@ using namespace vecgeom;
 
 int main(int argc, char* argv[]) {
   OPTION_INT(npoints,1024);
-  OPTION_INT(nrep,1024);
+  OPTION_INT(nrep,4);
   OPTION_DOUBLE(drmin,1.2);
   OPTION_DOUBLE(drmax,3.1);
   OPTION_DOUBLE(drtor,5.);

test/core/TestNavigationStatePool.cpp

@@ -27,7 +27,7 @@ int main()
   CudaManager::Instance().set_verbose(3);
   CudaManager::Instance().LoadGeometry();
 
-  // why to I have to do this??
+  // why do I have to do this here??
   CudaManager::Instance().Synchronize();
 
   CudaManager::Instance().PrintGeometry();

test/core/testVectorSafety.cpp

@@ -89,15 +89,14 @@ void testVectorNavigator( VPlacedVolume* world ){
    vecgeom::volumeUtilities::FillRandomDirections( dirs );
 
    // now setup all the navigation states
-   NavigationState ** states = new NavigationState*[np];
-   NavigationState ** newstates = new NavigationState*[np];
+   NavStatePool states(np, GeoManager::Instance().getMaxDepth() );
+   NavStatePool newstates(np, GeoManager::Instance().getMaxDepth() );
 
    vecgeom::SimpleNavigator nav;
    for (int i=0;i<np;++i){
-      // pSteps[i] = kInfinity;
-     pSteps[i] = (i%2)? 1 : VECGEOM_NAMESPACE::kInfinity;
-       states[i] =  NavigationState::MakeInstance( GeoManager::Instance().getMaxDepth() );
-       newstates[i] = NavigationState::MakeInstance( GeoManager::Instance().getMaxDepth() );
+       // pSteps[i] = kInfinity;
+       pSteps[i] = (i%2)? 1 : VECGEOM_NAMESPACE::kInfinity;
+       states[i]->Clear();
        nav.LocatePoint( world, points[i], *states[i], true);
    }
 

test/globalbenchmarks/LocatePointsBenchmark.cpp

@@ -45,7 +45,7 @@ void benchVecGeom( SOA3D<Precision> const & points, NavStatePool & statepool )
     Stopwatch timer;
     timer.Start();
     SimpleNavigator nav;
-    for(int i=0; i< points.size(); ++i){
+    for(unsigned int i=0; i< points.size(); ++i){
         nav.LocatePoint(
                 GeoManager::Instance().GetWorld(), points[i], *(statepool[i]), true );
     }
@@ -76,9 +76,9 @@ void benchGeant4( G4Navigator * nav, std::vector<G4ThreeVector> const & points,
     Stopwatch timer;
     timer.Start();
     G4VPhysicalVolume * vol;
-    for(int i=0; i < points.size(); ++i){
+    for(unsigned int i=0; i < points.size(); ++i){
         vol = nav->LocateGlobalPointAndSetup(points[i], NULL, false, true);
-
+        if(vol) ; // avoid compiler warning
         // this takes ages:
         // nav->LocateGlobalPointAndUpdateTouchable(points[i], pool[i], false);
     }
@@ -153,8 +153,8 @@ int main( int argc, char * argv[] )
     parser.Read( argv[2] );
 
     G4Navigator * nav = new G4Navigator();
-          nav->SetWorldVolume( const_cast<G4VPhysicalVolume *>(parser.GetWorldVolume()) );
-          nav->LocateGlobalPointAndSetup(G4ThreeVector(0,0,0), false );
+    nav->SetWorldVolume( const_cast<G4VPhysicalVolume *>(parser.GetWorldVolume()) );
+    nav->LocateGlobalPointAndSetup(G4ThreeVector(0,0,0), false );
     G4TouchableHistory ** Geant4statepool = new G4TouchableHistory*[npoints];
     for( int i=0;i<npoints;++i )
         {

test/globalbenchmarks/NavigationBenchmark.cpp

@@ -75,10 +75,10 @@ VPlacedVolume* SetupGeometry() {
 
 int main(int argc, char* argv[])
 {
-  OPTION_INT(npoints, 10000);
+  OPTION_INT(ntracks, 10000);
   OPTION_INT(nreps, 3);
   OPTION_STRING(geometry, "navBench");
-  OPTION_STRING(startVolName, "world");
+  OPTION_STRING(logvol, "world");
   OPTION_DOUBLE(bias, 0.8f);
 #ifdef VECGEOM_ROOT
   OPTION_BOOL(vis, false);
@@ -89,9 +89,8 @@ int main(int argc, char* argv[])
   OPTION_BOOL(help, false);
   if(help) return 0;
 
-  const VPlacedVolume *world = NULL;
   if(geometry.compare("navBench")==0) {
-    world = SetupGeometry();
+    const VPlacedVolume *world = SetupGeometry();
 
 #ifdef VECGEOM_ROOT
     // Exporting to ROOT file
@@ -118,7 +117,7 @@ int main(int argc, char* argv[])
   if(vis) {  // note that visualization block returns, excluding the rest of benchmark
     Visualizer visualizer;
     const VPlacedVolume* world = GeoManager::Instance().GetWorld();
-    world = GeoManager::Instance().FindPlacedVolume(startVolName.c_str());
+    world = GeoManager::Instance().FindPlacedVolume(logvol.c_str());
     visualizer.AddVolume( *world );
 
     Vector<Daughter> const* daughters = world->GetLogicalVolume()->daughtersp();
@@ -142,27 +141,26 @@ int main(int argc, char* argv[])
   }
 #endif
 
-  //testVectorSafety(world);
 
-  std::cout<<"\n*** Validating VecGeom navigation..."<< std::endl;
+  std::cout<<"\n*** Validating VecGeom navigation...\n";
 
-  const VPlacedVolume* startVolume = GeoManager::Instance().GetWorld();
-  if( startVolName.compare("world")!=0 ) {
-    startVolume = GeoManager::Instance().FindPlacedVolume(startVolName.c_str());
+  const LogicalVolume* startVolume = GeoManager::Instance().GetWorld()->GetLogicalVolume();
+  if( logvol.compare("world")!=0 ) {
+    startVolume = GeoManager::Instance().FindLogicalVolume(logvol.c_str());
   }
 
-  std::cout<<"NavigationBenchmark: startVolName=<"<< startVolName
-           <<">, startVolume="<< (startVolume ? startVolume->GetLabel() : NULL)
-           <<" - "<< *startVolume <<"\n";
+  std::cout<<"NavigationBenchmark: logvol=<"<< logvol
+           <<">, startVolume=<"<< (startVolume ? startVolume->GetLabel() : "NULL") <<">\n";
+  if(startVolume) std::cout<< *startVolume <<"\n";
+
+  int np = Min( ntracks, 1000 );  // no more than 1000 points used for validation
 
-  int np = Min( npoints, 1000 );  // no more than 1000 points used for validation
   SOA3D<Precision> points(np);
   SOA3D<Precision> dirs(np);
   SOA3D<Precision> locpts(np);
 
   vecgeom::volumeUtilities::FillGlobalPointsAndDirectionsForLogicalVolume(
-    startVolume->GetLogicalVolume(), locpts, points, dirs, bias, np);
-
+    startVolume, locpts, points, dirs, bias, np);
 
   // Must be validated before being benchmarked
   bool ok = validateVecGeomNavigation(np, points, dirs);
@@ -172,11 +170,11 @@ int main(int argc, char* argv[])
   }
   std::cout<<"VecGeom validation passed."<< std::endl;
 
-  // on mic.fnal.gov CPUs, loop execution takes ~70sec for npoints=10M
-  while(npoints<=10000) {
-    std::cout<<"\n*** Running navigation benchmarks with npoints="<<npoints<<" and nreps="<< nreps <<".\n";
-    runNavigationBenchmarks(startVolume, npoints, nreps, bias);
-    npoints*=10;
+  // on mic.fnal.gov CPUs, loop execution takes ~70sec for ntracks=10M
+  while(ntracks<=10000) {
+    std::cout<<"\n*** Running navigation benchmarks with ntracks="<<ntracks<<" and nreps="<< nreps <<".\n";
+    runNavigationBenchmarks(startVolume, ntracks, nreps, bias);
+    ntracks*=10;
   }
 
 

test/globalbenchmarks/NavigationBenchmarkFull.cpp

+/*
+ *  File: NavigationBenchmark.cpp
+ *
+ *  Created on: Oct 25, 2014
+ *      Author: swenzel, lima
+ */
+
+#ifdef VECGEOM_ROOT
+  #include "management/RootGeoManager.h"
+  #include "utilities/Visualizer.h"
+#endif
+
+#ifdef VECGEOM_GEANT4
+  #include "management/G4GeoManager.h"
+  #include "G4ThreeVector.hh"
+  // #include "G4TouchableHistoryHandle.hh"
+  #include "G4LogicalVolume.hh"
+  #include "G4LogicalVolumeStore.hh"
+  #include "G4PVPlacement.hh"
+  #include "G4GeometryManager.hh"
+#endif
+
+#include "benchmarking/NavigationBenchmarker.h"
+#include "test/benchmark/ArgParser.h"
+#include "volumes/utilities/VolumeUtilities.h"
+
+#include "management/GeoManager.h"
+#include "volumes/Box.h"
+#include "volumes/Orb.h"
+#include "volumes/Trapezoid.h"
+
+#include <fstream>
+
+using namespace VECGEOM_NAMESPACE;
+
+
+VPlacedVolume* SetupGeometry() {
+
+  UnplacedBox *worldUnplaced = new UnplacedBox(10, 10, 10);
+  UnplacedTrapezoid *trapUnplaced = new UnplacedTrapezoid(4,0,0,4,4,4,0,4,4,4,0);
+  UnplacedBox *boxUnplaced = new UnplacedBox(2.5, 2.5, 2.5);
+  UnplacedOrb *orbUnplaced = new UnplacedOrb(2.8);
+
+  LogicalVolume *world = new LogicalVolume("world",worldUnplaced);
+  LogicalVolume *trap  = new LogicalVolume("trap",trapUnplaced);
+  LogicalVolume *box   = new LogicalVolume("box",boxUnplaced);
+  LogicalVolume *orb   = new LogicalVolume("orb",orbUnplaced);
+
+  Transformation3D *ident = new Transformation3D( 0,  0,  0,  0,  0,  0);
+  orb->PlaceDaughter("orb1",box, ident);
+  trap->PlaceDaughter("box1",orb, ident);
+
+  Transformation3D *placement1 = new Transformation3D( 5,  5,  5,  0,  0,  0);
+  Transformation3D *placement2 = new Transformation3D(-5,  5,  5,  0,  0,  0); //45,  0,  0);
+  Transformation3D *placement3 = new Transformation3D( 5, -5,  5,  0,  0,  0); // 0, 45,  0);
+  Transformation3D *placement4 = new Transformation3D( 5,  5, -5,  0,  0,  0); // 0,  0, 45);
+  Transformation3D *placement5 = new Transformation3D(-5, -5,  5,  0,  0,  0); //45, 45,  0);
+  Transformation3D *placement6 = new Transformation3D(-5,  5, -5,  0,  0,  0); //45,  0, 45);
+  Transformation3D *placement7 = new Transformation3D( 5, -5, -5,  0,  0,  0); // 0, 45, 45);
+  Transformation3D *placement8 = new Transformation3D(-5, -5, -5,  0,  0,  0); //45, 45, 45);
+
+  world->PlaceDaughter("trap1",trap, placement1);
+  world->PlaceDaughter("trap2",trap, placement2);
+  world->PlaceDaughter("trap3",trap, placement3);
+  world->PlaceDaughter("trap4",trap, placement4);
+  world->PlaceDaughter("trap5",trap, placement5);
+  world->PlaceDaughter("trap6",trap, placement6);
+  world->PlaceDaughter("trap7",trap, placement7);
+  world->PlaceDaughter("trap8",trap, placement8);
+
+  VPlacedVolume* w = world->Place();
+  GeoManager::Instance().SetWorld(w);
+  GeoManager::Instance().CloseGeometry();
+  return w;
+}
+
+
+int main(int argc, char* argv[])
+{
+  OPTION_INT(ntracks, 10000);
+  // OPTION_INT(nreps, 3);
+  OPTION_STRING(geometry, "navBench");
+  OPTION_DOUBLE(bias, 0.0f);
+#ifdef VECGEOM_ROOT
+  OPTION_BOOL(vis, false);
+#endif
+
+  // default values used above are always printed.  If help true, stop now, so user will know which options
+  // are available, and what the default values are.
+  OPTION_BOOL(help, false);
+  if(help) return 0;
+
+  const VPlacedVolume *world = NULL;
+  if(world) ; // fix compilation warning
+  if(geometry.compare("navBench")==0) {
+    world = SetupGeometry();
+
+#ifdef VECGEOM_ROOT
+    // Exporting to ROOT file
+    RootGeoManager::Instance().ExportToROOTGeometry( world, "navBench.root" );
+    RootGeoManager::Instance().Clear();
+#endif
+  }
+
+  // Now try to read back in.  This is needed to make comparisons to VecGeom easily,
+  // since it builds VecGeom geometry based on the ROOT geometry and its TGeoNodes.
+#ifdef VECGEOM_ROOT
+  auto rootgeom = geometry+".root";
+  RootGeoManager::Instance().set_verbose(0);
+  RootGeoManager::Instance().LoadRootGeometry(rootgeom.c_str());
+#endif
+
+#ifdef VECGEOM_GEANT4
+  auto g4geom = geometry+".gdml";
+  G4GeoManager::Instance().LoadG4Geometry( g4geom.c_str() );
+#endif
+
+  std::string logvol("world");
+
+  // Visualization
+#ifdef VECGEOM_ROOT
+  if(vis) {  // note that visualization block returns, excluding the rest of benchmark
+    Visualizer visualizer;
+    const VPlacedVolume* world = GeoManager::Instance().GetWorld();
+    world = GeoManager::Instance().FindPlacedVolume(logvol.c_str());
+    visualizer.AddVolume( *world );
+
+    Vector<Daughter> const* daughters = world->GetLogicalVolume()->daughtersp();
+    for(int i=0; i<daughters->size(); ++i) {
+      VPlacedVolume const* daughter = (*daughters)[i];
+      Transformation3D const& trf1 = *(daughter->GetTransformation());
+      visualizer.AddVolume(*daughter, trf1);
+
+      // Vector<Daughter> const* daughters2 = daughter->GetLogicalVolume()->daughtersp();
+      // for(int ii=0; ii<daughters2->size(); ++ii) {
+      //   VPlacedVolume const* daughter2 = (*daughters2)[ii];
+      //   Transformation3D const& trf2 = *(daughter2->transformation());
+      //   Transformation3D comb = trf1;
+      //   comb.MultiplyFromRight(trf2);
+      //   visualizer.AddVolume(*daughter2, comb);
+      // }
+    }
+
+    visualizer.Show();
+    return 0;
+  }
+#endif
+
+  std::cout<<"\n*** Validating VecGeom navigation...\n";
+
+  const LogicalVolume* startVolume = GeoManager::Instance().GetWorld()->GetLogicalVolume();
+  //.. open input data file with volume names and loop over each volume
+  std::ifstream volumesFile("cms2015-logVolNames.lis");
+  if( !volumesFile.is_open() ) {
+    std::cout<<"\n*** Problems opening file with volume names!\n";
+    exit(-1);
+  }
+
+  std::string line;
+  while (getline(volumesFile,logvol)) {
+
+    std::cout<<"\n===================== volume: "<< logvol <<" =================\n";
+    // std::ofstream fout(logvol+".out");
+    // if(!fout.is_open()) {
+    //   std::cerr<<"Problems opening file: "<< logvol+".out\n";
+    //   continue;
+    // }
+
+    if( logvol.compare("world")!=0 ) {
+      startVolume = GeoManager::Instance().FindLogicalVolume(logvol.c_str());
+    }
+
+    std::cout<<"\nNavigationBenchmark: logvol=<"<< logvol
+             <<">, startVolume="<< (startVolume ? startVolume->GetLabel() : "NULL")
+             <<" - "<< *startVolume <<"\n";
+
+    int np = Min( ntracks, 1000 );  // no more than 1000 points used for validation
+    SOA3D<Precision> points(np);
+    SOA3D<Precision> dirs(np);
+    SOA3D<Precision> locpts(np);
+
+    vecgeom::volumeUtilities::FillGlobalPointsAndDirectionsForLogicalVolume(startVolume, locpts, points, dirs, bias, np);
+
+    // Must be validated before being benchmarked
+    bool ok = validateVecGeomNavigation(np, points, dirs);
+    if(!ok) {
+      std::cout<<"VecGeom validation failed."<< std::endl;
+      return 1;
+    }
+    std::cout<<"VecGeom validation passed."<< std::endl;
+
+    // on mic.fnal.gov CPUs, loop execution takes ~70sec for ntracks=10M
+    // while(ntracks<=1000) {
+    //    std::cout<<"\n*** Running navigation benchmarks with ntracks="<<ntracks<<" and nreps="<< nreps <<".\n";
+    //    runNavigationBenchmarks(startVolume, ntracks, nreps, bias);
+    //    ntracks*=10;
+    // }
+
+  }
+
+
+/*
+// GPU part
+  int nDevice;
+  cudaGetDeviceCount(&nDevice);
+
+  if(nDevice > 0) {
+    cudaDeviceReset();
+  }
+  else {
+    std::cout << "No Cuda Capable Device ... " << std::endl;
+    return 0;
+  }
+*/
+
+  return 0;
+}