R3BROOT
R3B analysis software
Loading...
Searching...
No Matches
R3BNeulandGeoCreator.cxx
Go to the documentation of this file.
2#include <FairGeoBuilder.h>
3#include <FairGeoInterface.h>
4#include <FairGeoLoader.h>
5#include <FairGeoMedia.h>
6#include <R3BException.h>
7#include <R3BLogger.h>
8#include <R3BNeulandCommon.h>
9#include <TGeoBBox.h>
10#include <TGeoCompositeShape.h>
11#include <TGeoCone.h>
12#include <TGeoManager.h>
13#include <TGeoMatrix.h>
14#include <TGeoShape.h>
15#include <TGeoVolume.h>
16#include <cmath>
17#include <cstdlib>
18#include <fmt/core.h>
19#include <string>
20#include <string_view>
21
23{
24 constexpr auto CONE_LENGHT_PADDING = 0.001; // cm
25 constexpr auto BAR_BASE_LENGTH = 125.0; // cm
26 constexpr auto BAR_CONE_LENGTH = 5.; // cm
27
28 const auto BC408_THICKNESS = 2.4; // cm
29 const auto BC408_CONE_RADIUS = 1.2; // cm
30 const auto AL_THICKNESS = 0.02; // cm
31 const auto TAPE_THICKNESS = 0.05; // cm
32
33 constexpr auto DEGREE_90 = 90.;
34 namespace
35 {
36 template <typename T, typename... Args>
37 auto TGeoManagerCreate(Args&&... args) -> T*
38 {
39 return std::make_unique<T>(std::forward<Args>(args)...).release();
40 }
41
42 auto create_cone_displacement(const std::string& name, double rotate_y_degree) -> TGeoCombiTrans*
43 {
44 auto rotation = TGeoRotation{};
45 rotation.RotateY(rotate_y_degree);
46 auto* transition = TGeoManagerCreate<TGeoCombiTrans>(
47 TGeoTranslation{ -(BAR_BASE_LENGTH + BAR_CONE_LENGTH), 0., 0. }, rotation);
48 transition->SetName(name.c_str());
49 transition->RegisterYourself();
50 return transition;
51 }
52
53 auto build_bar_shape(std::string_view name, const BarDimension& bar_dimension) -> TGeoShape*
54 {
55 auto* displace1 = create_cone_displacement("trc1", DEGREE_90);
56 auto* displace2 = create_cone_displacement("trc2", -DEGREE_90);
57
58 auto* main_bar_shape = TGeoManagerCreate<TGeoBBox>(
59 fmt::format("{}Box", name).c_str(), bar_dimension.length, bar_dimension.width, bar_dimension.width);
60 auto* cone_shape = TGeoManagerCreate<TGeoCone>(fmt::format("{}Cone", name).c_str(),
61 bar_dimension.cone_length + CONE_LENGHT_PADDING,
62 0.,
63 bar_dimension.cone_radius,
64 0.,
65 bar_dimension.width * std::sqrt(2.));
66 auto* bar_tail_shape = TGeoManagerCreate<TGeoBBox>(fmt::format("{}ConeBox", name).c_str(),
67 bar_dimension.width,
68 bar_dimension.width,
69 bar_dimension.cone_length);
70 return TGeoManagerCreate<TGeoCompositeShape>(name.data(),
71 fmt::format("{}+(({}*{}):{})+(({}*{}):{})",
72 main_bar_shape->GetName(),
73 bar_tail_shape->GetName(),
74 cone_shape->GetName(),
75 displace1->GetName(),
76 bar_tail_shape->GetName(),
77 cone_shape->GetName(),
78 displace2->GetName())
79 .c_str());
80 }
81
82 } // namespace
83
84 auto Creator::construct_volume(int num_of_planes, FairGeoLoader* geo_loader) -> TGeoVolume*
85 {
86 read_material_from_file(geo_loader);
90 return build_detector(num_of_planes);
91 }
92
93 auto Creator::build_detector(int num_of_planes) -> TGeoVolume*
94 {
95 auto* neuland = TGeoManagerCreate<TGeoVolumeAssembly>("volNeuland");
96
97 auto total_module_size = num_of_planes * BarsPerPlane;
98 const auto first_plane_front_z = -num_of_planes * BarSize_Z / 2.;
99
100 auto rot_zero = TGeoRotation{};
101 auto rot_90 = TGeoRotation{};
102 rot_90.RotateZ(DEGREE_90);
103
104 for (int module_id{}; module_id < total_module_size; ++module_id)
105 {
106 const auto plane_id = ModuleID2PlaneID(module_id);
107 const auto is_horizontal = IsPlaneIDHorizontal(plane_id);
108 const auto z_pos = PlaneID2ZPos(plane_id) + first_plane_front_z;
109 const auto vertical_displacement = GetBarVerticalDisplacement(module_id + 1);
110 auto* translation_rotation = TGeoManagerCreate<TGeoCombiTrans>();
111 if (is_horizontal)
112 {
113 translation_rotation->SetTranslation(0, vertical_displacement, z_pos);
114 translation_rotation->SetRotation(rot_zero);
115 }
116 else
117 {
118 translation_rotation->SetTranslation(vertical_displacement, 0, z_pos);
119 translation_rotation->SetRotation(rot_90);
120 }
121 neuland->AddNode(bar_, module_id + 1, translation_rotation);
122 }
123
124 return neuland;
125 }
126
134
136 {
138 material_poly_ = build_material("polyethylene");
139 material_Al_ = build_material("aluminium");
140 }
141
143 {
144 auto dimension = BarDimension{};
145
146 dimension.length = BAR_BASE_LENGTH;
147 dimension.width = BC408_THICKNESS;
148 dimension.cone_radius = BC408_CONE_RADIUS;
149 dimension.cone_length = BAR_CONE_LENGTH;
150 shape_scintillator_ = build_bar_shape("shapeBC408", dimension);
151
152 dimension.width += AL_THICKNESS;
153 dimension.cone_radius += AL_THICKNESS;
154 auto* shape_Al_solid = build_bar_shape("shapeAlWrappingSolid", dimension);
155 shape_Al_wrapping_ = TGeoManagerCreate<TGeoCompositeShape>(
156 "shapeAlWrapping",
157 fmt::format("{} - {}", shape_Al_solid->GetName(), shape_scintillator_->GetName()).c_str());
158
159 dimension.width += TAPE_THICKNESS;
160 dimension.cone_radius += TAPE_THICKNESS;
161 auto* shape_tape_solid = build_bar_shape("shapeTapeWrappingSolid", dimension);
162 shape_tape_wrapping_ = TGeoManagerCreate<TGeoCompositeShape>(
163 "shapeTapeWrapping",
164 fmt::format("{} - {}", shape_tape_solid->GetName(), shape_Al_solid->GetName()).c_str());
165 }
166
167 auto Creator::build_bar_volume() -> TGeoVolume*
168 {
169 auto* bar = TGeoManagerCreate<TGeoVolumeAssembly>("volPaddle");
170 bar->AddNode(scintillator_, 1);
171 bar->AddNode(Al_wrapping_, 1);
172 bar->AddNode(tape_wrapping_, 1);
173 return bar;
174 }
175
176 auto Creator::build_scintillator() -> TGeoVolume*
177 {
178 auto* scintillator = TGeoManagerCreate<TGeoVolume>("volBC408", shape_scintillator_, material_BC408_);
179 static constexpr auto greyish_blue = 33;
180 static constexpr auto transparency = 30;
181 scintillator->SetLineColor(greyish_blue);
182 scintillator->SetTransparency(transparency);
183 return scintillator;
184 }
185
186 auto Creator::build_Al_wrapping() -> TGeoVolume*
187 {
188 auto* Al_wrapping = TGeoManagerCreate<TGeoVolume>("volAlWrapping", shape_Al_wrapping_, material_Al_);
189 static constexpr auto greyish_silver = 17;
190 Al_wrapping->SetLineColor(greyish_silver); // grey/silver
191 return Al_wrapping;
192 }
193
194 auto Creator::build_tape_wrapping() -> TGeoVolume*
195 {
196 auto* tape_wrapping = TGeoManagerCreate<TGeoVolume>("volTapeWrapping", shape_tape_wrapping_, material_poly_);
197 static constexpr auto black = 1;
198 tape_wrapping->SetLineColor(black); // grey/silver
199 return tape_wrapping;
200 }
201
202 auto Creator::build_material(const std::string& material_name) -> TGeoMedium*
203 {
204 auto* fair_medium = geo_media_->getMedium(material_name.c_str());
205 if (fair_medium == nullptr)
206 {
207 R3BLOG(error, fmt::format("FairGeoMedium {} not found!", material_name));
208 }
209
210 if (geo_builder_ == nullptr)
211 {
212 throw R3B::runtime_error("geo_builder is nullptr!");
213 }
214 geo_builder_->createMedium(fair_medium);
215 if (gGeoManager == nullptr)
216 {
217 throw R3B::runtime_error("gGeoManager is nullptr!");
218 }
219 auto* material = gGeoManager->GetMedium(material_name.c_str());
220 if (material == nullptr)
221 {
222 R3BLOG(error, fmt::format("TGeoMedium {} not found!", material_name));
223 }
224 return material;
225 }
226
227 void Creator::read_material_from_file(FairGeoLoader* geo_loader)
228 {
229 auto* geo_interface = geo_loader->getGeoInterface();
230 const auto* working_dir = std::getenv("VMCWORKDIR");
231 if (working_dir == nullptr)
232 {
233 throw R3B::logic_error("Environment variable \"VMCWORKDIR\" is not defined!");
234 }
235 geo_interface->setMediaFile(fmt::format("{}/geometry/media_r3b.geo", working_dir).c_str());
236 geo_interface->readMedia();
237
238 geo_builder_ = geo_loader->getGeoBuilder();
239 geo_media_ = geo_interface->getMedia();
240 }
241} // namespace R3B::Neuland::Geometry
#define R3BLOG(severity, x)
Definition R3BLogger.h:32
auto construct_volume(int num_of_planes, FairGeoLoader *geo_loader) -> TGeoVolume *
auto build_detector(int num_of_planes) -> TGeoVolume *
void read_material_from_file(FairGeoLoader *geo_loader)
auto build_material(const std::string &material) -> TGeoMedium *
constexpr auto BarsPerPlane
constexpr auto ModuleID2PlaneID(int moduleID) -> int
constexpr auto GetBarVerticalDisplacement(int module_num) -> double
constexpr auto BarSize_Z
constexpr auto IsPlaneIDHorizontal(int plane_id) -> bool
constexpr auto PlaneID2ZPos(int plane_id) -> T