File: GenerateSurfacePoints.cpp

package info (click to toggle)
vecgeom 1.2.8%2Bdfsg-2
  • links: PTS, VCS
  • area: main
  • in suites: sid, trixie
  • size: 24,016 kB
  • sloc: cpp: 88,803; ansic: 6,888; python: 1,035; sh: 582; sql: 538; makefile: 23
file content (92 lines) | stat: -rw-r--r-- 2,987 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
#include "VecGeomTest/RootGeoManager.h"

#include "VecGeom/volumes/LogicalVolume.h"
#include "VecGeom/volumes/PlacedVolume.h"
#include "VecGeom/volumes/UnplacedBox.h"
#include "VecGeomTest/Benchmarker.h"
#include "VecGeom/management/GeoManager.h"
#include "VecGeom/volumes/UnplacedBox.h"
#include "TGeoManager.h"
#include "TGeoVolume.h"
#include "TGeoBBox.h"
#include "VecGeomTest/Visualizer.h"
#include <string>
#include <cmath>
#include <iostream>

#ifdef VECGEOM_GEANT4
#include "G4ThreeVector.hh"
#include "G4VSolid.hh"
#endif

size_t N = 1000;

using namespace vecgeom;

int main(int argc, char *argv[])
{

  if (argc < 3) {
    std::cerr << "need to give root geometry file + logical volume name\n";
    std::cerr << "example: " << argv[0] << " cms2015.root CALO\n";
    return 1;
  }

  TGeoManager::Import(argv[1]);
  std::string testvolume(argv[2]);
  int found               = 0;
  TGeoVolume *foundvolume = gGeoManager->FindVolumeFast(argv[2]);
  if (!foundvolume) {
    // now try to find shape with logical volume name given on the command line
    TObjArray *vlist = gGeoManager->GetListOfVolumes();
    for (auto i = 0; i < vlist->GetEntries(); ++i) {
      TGeoVolume *vol = reinterpret_cast<TGeoVolume *>(vlist->At(i));
      std::string fullname(vol->GetName());

      std::size_t founds = fullname.compare(testvolume);
      if (founds == 0) {
        found++;
        foundvolume = vol;
        std::cerr << "found matching volume " << fullname << " of type " << vol->GetShape()->ClassName() << "\n";
      }
    }
    std::cerr << "volume found " << found << " times \n";
    foundvolume->GetShape()->InspectShape();
    std::cerr << "volume capacity " << foundvolume->GetShape()->Capacity() << "\n";
  }

  // now get the VecGeom shape and benchmark it
  if (foundvolume) {
    LogicalVolume *converted     = RootGeoManager::Instance().Convert(foundvolume);
    VPlacedVolume *vecgeomplaced = converted->Place();

    // generate POINTS ON SURFACE FROM VECGEOM
    for (size_t i = 0; i < N; ++i) {
      auto sp = vecgeomplaced->GetUnplacedVolume()->SamplePointOnSurface();
      std::cerr << "SPVG " << sp.x() << " " << sp.y() << " " << sp.z() << "\n";
    }

// generate POINTS ON SURFACE FOR G4
#ifdef VECGEOM_GEANT4
    auto g4volume = vecgeomplaced->ConvertToGeant4();
    for (size_t i = 0; i < N; ++i) {
      auto sp = g4volume->GetPointOnSurface();
      std::cerr << "SPG4 " << sp.x() << " " << sp.y() << " " << sp.z() << "\n";
    }

#endif
    // generate POINTS ON SURFACE FOR ROOT
    double points[3 * N];
    bool success = foundvolume->GetShape()->GetPointsOnSegments(N, points);
    if (success) {
      for (size_t i = 0; i < N; ++i) {
        std::cerr << "SPR " << points[3 * i] << " " << points[3 * i + 1] << " " << points[3 * i + 2] << "\n";
      }
    } else {
      std::cerr << "ROOT cannot generate surface points \n";
    }
  } else {
    std::cerr << " NO SUCH VOLUME [" << testvolume << "] FOUND ... EXITING \n";
    return 1;
  }
}