LDMX Software
EcalTrackingGeometry.cxx
1#include "Tracking/geo/EcalTrackingGeometry.h"
2
3namespace tracking::geo {
4
5EcalTrackingGeometry::EcalTrackingGeometry(std::string gdmlfile,
6 Acts::GeometryContext* gctx,
7 bool debug) {
8 _debug = debug;
9 gctx_ = gctx;
10
11 // Build The rotation matrix to the tracking frame
12 // Rotate the sensors to be orthogonal to X
13 double rotationAngle = M_PI * 0.5;
14
15 // 0 0 -1
16 // 0 1 0
17 // 1 0 0
18
19 // This rotation is needed to have the plane orthogonal to the X direction.
20 // Rotation of the surfaces
21 Acts::Vector3 xPos1(cos(rotationAngle), 0., sin(rotationAngle));
22 Acts::Vector3 yPos1(0., 1., 0.);
23 Acts::Vector3 zPos1(-sin(rotationAngle), 0., cos(rotationAngle));
24
25 y_rot_.col(0) = xPos1;
26 y_rot_.col(1) = yPos1;
27 y_rot_.col(2) = zPos1;
28
29 // Rotate the sensors to put them in the proper orientation in Z
30 Acts::Vector3 xPos2(1., 0., 0.);
31 Acts::Vector3 yPos2(0., cos(rotationAngle), sin(rotationAngle));
32 Acts::Vector3 zPos2(0., -sin(rotationAngle), cos(rotationAngle));
33
34 x_rot_.col(0) = xPos2;
35 x_rot_.col(1) = yPos2;
36 x_rot_.col(2) = zPos2;
37
38 G4GDMLParser parser;
39
40 // Validation requires internet
41 parser.Read(
42 "/Users/pbutti/sw/ldmx-sw/Detectors/data/ldmx-det-v12/detector.gdml",
43 false);
44 G4VPhysicalVolume* fWorldPhysVol = parser.GetWorldVolume();
45 if (debug) {
46 std::cout << "World position" << std::endl;
47 std::cout << fWorldPhysVol->GetTranslation() << std::endl;
48 }
49
50 if (debug) std::cout << "Loop on world daughters" << std::endl;
51
52 _ECAL = findDaughterByName(fWorldPhysVol, "em_calorimeters_PV");
53
54 if (_ECAL) {
55 std::cout << _ECAL->GetName() << std::endl;
56 std::cout << _ECAL->GetTranslation() << std::endl;
57 std::cout << _ECAL->GetLogicalVolume()->GetNoDaughters() << std::endl;
58 }
59
60 // getAllDaughters(_ECAL);
61
62 if (_debug) {
63 std::cout << "Looking for later components" << std::endl;
64 }
65
66 // Layers Config
67 std::vector<Acts::CuboidVolumeBuilder::LayerConfig> layersCfg;
68
69 // Get the layers
70 std::vector<std::reference_wrapper<G4VPhysicalVolume>> components;
71 std::vector<std::shared_ptr<const Acts::Surface>> sensor_surfaces;
72
73 // PS
74 getComponentRing("_a_", "Si", components);
75
76 for (auto& component : components) {
77 sensor_surfaces.push_back(convertHexToActsSurface(component));
78 }
79
80 for (int i = 0; i < sensor_surfaces.size(); i += 14) {
81 std::vector<std::shared_ptr<const Acts::Surface>> rings;
82 for (int j = 0; j < 14; j++) {
83 rings.push_back(sensor_surfaces.at(i + j));
84 }
85
86 Acts::CuboidVolumeBuilder::LayerConfig lyCfg = buildLayerConfig(rings);
87 layersCfg.push_back(lyCfg);
88 }
89
90 components.clear();
91 sensor_surfaces.clear();
92
93 // A
94
95 getComponentRing("_a_", "Si", components);
96
97 for (auto& component : components) {
98 sensor_surfaces.push_back(convertHexToActsSurface(component));
99 }
100
101 for (int i = 0; i < sensor_surfaces.size(); i += 14) {
102 std::vector<std::shared_ptr<const Acts::Surface>> rings;
103 for (int j = 0; j < 14; j++) {
104 rings.push_back(sensor_surfaces.at(i + j));
105 }
106
107 Acts::CuboidVolumeBuilder::LayerConfig lyCfg = buildLayerConfig(rings);
108 layersCfg.push_back(lyCfg);
109 }
110
111 components.clear();
112 sensor_surfaces.clear();
113
114 /*
115 //B
116
117
118 getComponentRing("_b_", "Si",components);
119
120 for (auto& component : components) {
121 sensor_surfaces.push_back(convertHexToActsSurface(component));
122 }
123
124 for (int i = 0; i< sensor_surfaces.size(); i+=14) {
125
126 std::vector<std::shared_ptr<const Acts::Surface> > rings;
127 for (int j = 0; j < 14 ; j++) {
128 rings.push_back(sensor_surfaces.at(i+j));
129 }
130
131 Acts::CuboidVolumeBuilder::LayerConfig lyCfg = buildLayerConfig(rings);
132 layersCfg.push_back(lyCfg);
133 }
134
135 components.clear();
136 sensor_surfaces.clear();
137
138
139
140 //C
141
142 std::cout<<"Getting the components for C"<<std::endl;
143 getComponentRing("_c_", "Si",components);
144
145 std::cout<<"Building the surfaces"<<std::endl;
146 for (auto& component : components) {
147 sensor_surfaces.push_back(convertHexToActsSurface(component));
148 }
149
150 std::cout<<"Selecting the rings"<<std::endl;
151
152 for (int i = 0; i< sensor_surfaces.size(); i+=14) {
153
154 std::vector<std::shared_ptr<const Acts::Surface> > rings;
155
156 std::cout<<"Getting ring " << i / 14 << std::endl;
157 for (int j = 0; j < 14 ; j++) {
158 rings.push_back(sensor_surfaces.at(i+j));
159 }
160
161 Acts::CuboidVolumeBuilder::LayerConfig lyCfg = buildLayerConfig(rings);
162 layersCfg.push_back(lyCfg);
163 }
164
165 components.clear();
166 sensor_surfaces.clear();
167
168
169 //D
170
171 getComponentRing("_d_", "Si",components);
172
173 for (auto& component : components) {
174 sensor_surfaces.push_back(convertHexToActsSurface(component));
175 }
176
177 for (int i = 0; i< sensor_surfaces.size(); i+=14) {
178
179 std::vector<std::shared_ptr<const Acts::Surface> > rings;
180 for (int j = 0; j < 14 ; j++) {
181 rings.push_back(sensor_surfaces.at(i+j));
182 }
183
184 Acts::CuboidVolumeBuilder::LayerConfig lyCfg = buildLayerConfig(rings);
185 layersCfg.push_back(lyCfg);
186 }
187
188 components.clear();
189 sensor_surfaces.clear();
190 */
191
192 Acts::CuboidVolumeBuilder::VolumeConfig ecal_vcfg;
193
194 // TODO change this with the translation!
195 ecal_vcfg.position = {0., 0., 0.};
196 G4Box* ecal_box = (G4Box*)_ECAL->GetLogicalVolume()->GetSolid();
197 // ecal_vcfg.length = {ecal_box->GetXHalfLength() * 2.,
198 // ecal_box->GetYHalfLength() * 2. , ecal_box->GetZHalfLength() * 2.};
199
200 ecal_vcfg.length = {ecal_box->GetZHalfLength() * 2.,
201 ecal_box->GetXHalfLength() * 2.,
202 ecal_box->GetYHalfLength() * 2.};
203
204 ecal_vcfg.name = "Ecal_volume";
205 ecal_vcfg.layerCfg = layersCfg;
206
207 // Initialize the Cuboid Volume Builder
208 Acts::CuboidVolumeBuilder cvb;
209 Acts::CuboidVolumeBuilder::Config config;
210 config.position = {0., 0., 0.};
211 config.length = {2000, 2000, 2000};
212 config.volumeCfg = {ecal_vcfg};
213
214 cvb.setConfig(config);
215
216 Acts::TrackingGeometryBuilder::Config tgbCfg;
217 tgbCfg.trackingVolumeBuilders.push_back(
218 [=](const auto& cxt, const auto& inner, const auto&) {
219 return cvb.trackingVolume(cxt, inner, nullptr);
220 });
221 Acts::TrackingGeometryBuilder tgb(tgbCfg);
222 tGeometry_ = tgb.trackingGeometry(*gctx_);
223
224 dumpGeometry("./ecal_test/");
225}
226
227G4VPhysicalVolume* EcalTrackingGeometry::findDaughterByName(
228 G4VPhysicalVolume* pvol, G4String name) {
229 G4LogicalVolume* lvol = pvol->GetLogicalVolume();
230 for (G4int i = 0; i < lvol->GetNoDaughters(); i++) {
231 G4VPhysicalVolume* fDaughterPhysVol = lvol->GetDaughter(i);
232 if (fDaughterPhysVol->GetName() == name) return fDaughterPhysVol;
233 }
234
235 return nullptr;
236}
237
238// A super layer is a full layer composed by multiple duouble-silicon layers.
239// One pass a super layer, and this function returns a vector of all the
240// components of the super layer of a certain type It then orders then by z
241// location And returns a map where at every z the ring is organised in a vector
242// with the numbering scheme as described here:
243// https://www.overleaf.com/project/5f3d5a7e36c7f90001eddfc6
244
245void EcalTrackingGeometry::getComponentRing(
246 std::string super_layer_name, std::string component_type,
247 std::vector<std::reference_wrapper<G4VPhysicalVolume>>& components) {
248 G4LogicalVolume* l_ECAL = _ECAL->GetLogicalVolume();
249 for (G4int i = 0; i < l_ECAL->GetNoDaughters(); i++) {
250 // std::cout<<"DEBUG:: Getting daughter.."<< i << " of "<<
251 // l_ECAL->GetNoDaughters()<<std::endl;
252 std::string sln = l_ECAL->GetDaughter(i)->GetName();
253
254 // std::cout<<sln<<std::endl;
255
256 // Look for the specific super layer of the specific component
257 if (sln.find(super_layer_name) != std::string::npos &&
258 sln.find(component_type) != std::string::npos) {
259 // std::cout<<"Found and adding " << sln<<std::endl;
260 components.push_back(*(l_ECAL->GetDaughter(i)));
261 }
262 }
263
264 // Sort the components along z
265 // std::cout<<"Sorting..." <<std::endl;
266 sort(components.begin(), components.end(), compareZlocation);
267}
268
269void EcalTrackingGeometry::getAllDaughters(G4VPhysicalVolume* pvol) {
270 G4LogicalVolume* lvol = pvol->GetLogicalVolume();
271
272 for (G4int i = 0; i < lvol->GetNoDaughters(); i++) {
273 G4VPhysicalVolume* fDaughterPhysVol = lvol->GetDaughter(i);
274 if (_debug) {
275 std::cout << "name::" << fDaughterPhysVol->GetName() << std::endl;
276 std::cout << "pos::" << fDaughterPhysVol->GetTranslation() << std::endl;
277 std::cout << "n_dau::"
278 << fDaughterPhysVol->GetLogicalVolume()->GetNoDaughters()
279 << std::endl;
280 std::cout << "replica::" << fDaughterPhysVol->IsReplicated() << std::endl;
281 std::cout << "copyNR::" << fDaughterPhysVol->GetCopyNo() << std::endl;
282 }
283 }
284}
285
286void EcalTrackingGeometry::dumpGeometry(const std::string& outputDir) {
287 // Should fail if already exists
288 boost::filesystem::create_directory(outputDir);
289
290 double outputScalor = 1.0;
291 size_t outputPrecision = 6;
292
293 Acts::ObjVisualization3D objVis(outputPrecision, outputScalor);
294 Acts::ViewConfig containerView = Acts::ViewConfig({220, 220, 220});
295 Acts::ViewConfig volumeView = Acts::ViewConfig({220, 220, 0});
296 Acts::ViewConfig sensitiveView = Acts::ViewConfig({0, 180, 240});
297 Acts::ViewConfig passiveView = Acts::ViewConfig({240, 280, 0});
298 Acts::ViewConfig gridView = Acts::ViewConfig({220, 0, 0});
299
300 Acts::GeometryView3D::drawTrackingVolume(
301 objVis, *(tGeometry_->highestTrackingVolume()), *gctx_, containerView,
302 volumeView, passiveView, sensitiveView, gridView, true, "", ".");
303}
304
305std::shared_ptr<const Acts::Surface>
306EcalTrackingGeometry::convertHexToActsSurface(const G4VPhysicalVolume& phex) {
307 G4Polyhedra* hex = (G4Polyhedra*)phex.GetLogicalVolume()->GetSolid();
308 G4Material* hex_mat = phex.GetLogicalVolume()->GetMaterial();
309
310 double side = hex->GetCorner(3).r;
311 double height = (side / 2.) * sqrt(3);
312 double thickness = 0.5 * Acts::UnitConstants::mm;
313
314 std::shared_ptr<Acts::PlaneSurface> surface;
315
323 std::shared_ptr<const Acts::DiamondBounds> d_bounds =
324 std::make_shared<Acts::DiamondBounds>(
325 side / 2. * Acts::UnitConstants::mm, side * Acts::UnitConstants::mm,
326 side / 2. * Acts::UnitConstants::mm, height * Acts::UnitConstants::mm,
327 height * Acts::UnitConstants::mm);
328
329 // Form the position and rotation
330 // Acts::Vector3 pos(phex->GetTranslation().x(), phex->GetTranslation().y(),
331 // phex->GetTranslation().z());
332
333 Acts::Vector3 pos(phex.GetTranslation().z(), phex.GetTranslation().x(),
334 phex.GetTranslation().y());
335
336 Acts::RotationMatrix3 rotation;
337
338 // ConvertG4Rot(phex->GetObjectRotation(),rotation);
339 ConvertG4Rot(phex.GetObjectRotationValue(), rotation);
340
341 // Now rotate to the tracking frame
342 rotation = x_rot_ * y_rot_ * rotation;
343
344 Acts::Translation3 translation(pos);
345 Acts::Transform3 transform(translation * rotation);
346 surface = Acts::Surface::makeShared<Acts::PlaneSurface>(transform, d_bounds);
347
348 Acts::Material silicon = Acts::Material::fromMassDensity(
349 hex_mat->GetRadlen(), hex_mat->GetNuclearInterLength(), hex_mat->GetA(),
350 hex_mat->GetZ(), hex_mat->GetDensity());
351
352 Acts::MaterialSlab silicon_slab(silicon, thickness);
353 surface->assignSurfaceMaterial(
354 std::make_shared<Acts::HomogeneousSurfaceMaterial>(silicon_slab));
355
356 // surface->toStreamImpl(*gctx_, std::cout);
357 surface->toStream(*gctx_);
358
359 return surface;
360}
361
362// Convert rotation
363
364void EcalTrackingGeometry::ConvertG4Rot(const G4RotationMatrix& g4rot,
365 Acts::RotationMatrix3& rot) {
366 rot(0, 0) = g4rot.xx();
367 rot(0, 1) = g4rot.xy();
368 rot(0, 2) = g4rot.xz();
369
370 rot(1, 0) = g4rot.yx();
371 rot(1, 1) = g4rot.yy();
372 rot(1, 2) = g4rot.yz();
373
374 rot(2, 0) = g4rot.zx();
375 rot(2, 1) = g4rot.zy();
376 rot(2, 2) = g4rot.zz();
377}
378
379// Convert translation
380
381Acts::Vector3 EcalTrackingGeometry::ConvertG4Pos(const G4ThreeVector& g4pos) {
382 std::cout << std::endl;
383 std::cout << "g4pos::" << g4pos << std::endl;
384 Acts::Vector3 trans{g4pos.x(), g4pos.y(), g4pos.z()};
385 std::cout << trans << std::endl;
386 return trans;
387}
388
389// A layer is formed by 2 sides of silicon + material
390//
391Acts::CuboidVolumeBuilder::LayerConfig EcalTrackingGeometry::buildLayerConfig(
392 const std::vector<std::shared_ptr<const Acts::Surface>>& rings,
393 double clearance, bool active) {
394 Acts::CuboidVolumeBuilder::LayerConfig lcfg;
395
396 for (auto& hex : rings) {
397 lcfg.surfaces.push_back(hex);
398 }
399
400 // Get the surface thickness
401 Acts::Vector2 ploc{0., 0.};
402 double thickness =
403 rings.front()->surfaceMaterial()->materialSlab(ploc).thickness();
404
405 lcfg.envelopeX = std::array<double, 2>{thickness / 2. + clearance,
406 thickness / 2. + clearance};
407 lcfg.active = active;
408
409 return lcfg;
410}
411} // namespace tracking::geo
std::shared_ptr< const Acts::Surface > convertHexToActsSurface(const G4VPhysicalVolume &phex)
Visualization.