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
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
|
/*
* BenchmarkShapeFromROOTFile.cpp
*
* Created on: Jan 26, 2015
* Author: swenzel
*/
#include "VecGeomTest/RootGeoManager.h"
#include "VecGeom/volumes/LogicalVolume.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>
using namespace vecgeom;
// benchmarking any available shape (logical volume) found in a ROOT file
// usage: BenchmarkShapeFromROOTFile detector.root logicalvolumename
// logicalvolumename should not contain trailing pointer information
int main(int argc, char *argv[])
{
if (argc < 3) {
std::cerr << "need to give root geometry file and logical volume name";
}
TGeoManager::Import(argv[1]);
std::string testvolume(argv[2]);
int found = 0;
TGeoVolume *foundvolume = NULL;
// 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);
converted->Print();
// get the bounding box
TGeoBBox *bbox = dynamic_cast<TGeoBBox *>(foundvolume->GetShape());
double bx, by, bz;
double const *offset;
if (bbox) {
bx = bbox->GetDX();
by = bbox->GetDY();
bz = bbox->GetDZ();
bbox->InspectShape();
offset = bbox->GetOrigin();
UnplacedBox worldUnplaced((bx + fabs(offset[0])) * 4, (by + fabs(offset[1])) * 4, (bz + fabs(offset[2])) * 4);
LogicalVolume world("world", &worldUnplaced);
// for the moment at origin
Transformation3D placement(0, 0, 0);
VPlacedVolume const *vol = world.PlaceDaughter("testshape", converted, &placement);
VPlacedVolume *worldPlaced = world.Place();
GeoManager::Instance().SetWorld(worldPlaced);
Benchmarker tester(GeoManager::Instance().GetWorld());
tester.SetTolerance(1E-6);
tester.SetVerbosity(3);
tester.SetPoolMultiplier(1);
tester.SetRepetitions(1);
tester.SetPointCount(100);
tester.CompareMetaInformation();
int returncodeIns = tester.RunInsideBenchmark();
if (returncodeIns) {
// try to visualize if there was a problem
Visualizer visualizer;
visualizer.AddVolume(*vol);
if (tester.GetProblematicContainPoints().size() > 0) {
for (auto v : tester.GetProblematicContainPoints()) {
visualizer.AddPoint(v);
// for debugging purpose
std::cerr << " " << vol->Contains(v) << "\n";
std::cout << v << "\n";
}
visualizer.Show();
}
}
int returncodeToIn = tester.RunToInBenchmark();
if (returncodeToIn) {
Visualizer visualizer;
visualizer.AddVolume(*vol);
for (auto v : tester.GetProblematicRays()) {
auto point = v.first; // first and second are direct member of std::pair ...
visualizer.AddPoint(point);
visualizer.AddLine(v.first, v.first + std::max(by, std::max(bz, bx)) * v.second);
}
visualizer.Show();
}
int returncodeToOut = tester.RunToOutBenchmark();
if (returncodeToOut) {
Visualizer visualizer;
visualizer.AddVolume(*vol);
for (auto v : tester.GetProblematicRays()) {
auto point = v.first; // first and second are direct member of std::pair ...
visualizer.AddPoint(point);
visualizer.AddLine(v.first, v.first + std::sqrt(4 * (worldUnplaced.x() * worldUnplaced.x() +
worldUnplaced.y() * worldUnplaced.y() +
worldUnplaced.z() * worldUnplaced.z())));
}
visualizer.Show();
}
return returncodeIns + returncodeToIn + returncodeToOut;
}
return 1;
} else {
std::cerr << " NO SUCH VOLUME [" << testvolume << "] FOUND ... EXITING \n";
return 1;
}
}
|
