File indexing completed on 2024-06-18 07:05:20
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014 #include "DD4hep/DetFactoryHelper.h"
0015 #include "XML/Layering.h"
0016 #include "Math/Point2D.h"
0017 #include "TGeoPolygon.h"
0018
0019 using namespace std;
0020 using namespace dd4hep;
0021 using namespace dd4hep::detail;
0022
0023 typedef ROOT::Math::XYPoint Point;
0024
0025 vector<vector<Point>> fiberPositions(double radius, double x_spacing, double z_spacing,
0026 double x, double z, double phi, double spacing_tol = 1e-2);
0027 std::pair<int, int> getNdivisions(double x, double z, double dx, double dz);
0028 vector<tuple<int, Point, Point, Point, Point>> gridPoints(int div_x, int div_z, double x, double z, double phi);
0029
0030
0031 void buildFibers(Detector& desc, SensitiveDetector &sens, Volume &mother, xml_comp_t x_fiber,
0032 const std::tuple<double, double, double, double> &dimensions);
0033 void buildSupport(Detector& desc, Volume &mother, xml_comp_t x_support,
0034 const std::tuple<double, double, double, double> &dimensions);
0035
0036
0037
0038 static Ref_t create_detector(Detector& desc, xml_h e, SensitiveDetector sens) {
0039 Layering layering (e);
0040 xml_det_t x_det = e;
0041 Material air = desc.air();
0042 int det_id = x_det.id();
0043 string det_name = x_det.nameStr();
0044 double offset = x_det.attr<double>(_Unicode(offset));
0045 xml_comp_t x_dim = x_det.dimensions();
0046 int nsides = x_dim.numsides();
0047 double inner_r = x_dim.rmin();
0048 double dphi = (2*M_PI/nsides);
0049 double hphi = dphi/2;
0050
0051 DetElement sdet (det_name, det_id);
0052 Volume motherVol = desc.pickMotherVolume(sdet);
0053
0054 Assembly envelope (det_name);
0055 Transform3D tr = Translation3D(0, 0, offset) * RotationZ(hphi);
0056 PlacedVolume env_phv = motherVol.placeVolume(envelope, tr);
0057 sens.setType("calorimeter");
0058
0059 env_phv.addPhysVolID("system",det_id);
0060 sdet.setPlacement(env_phv);
0061
0062
0063 DetElement stave_det("stave0", det_id);
0064 Assembly mod_vol("stave");
0065
0066
0067 double l_pos_z = inner_r;
0068 {
0069
0070 double tan_hphi = std::tan(hphi);
0071 double l_dim_y = x_dim.z()/2.;
0072
0073
0074 int l_num = 1;
0075 for(xml_coll_t li(x_det, _U(layer)); li; ++li) {
0076 xml_comp_t x_layer = li;
0077 int repeat = x_layer.repeat();
0078 double l_space_between = dd4hep::getAttrOrDefault(x_layer, _Unicode(space_between), 0.);
0079 double l_space_before = dd4hep::getAttrOrDefault(x_layer, _Unicode(space_before), 0.);
0080 l_pos_z += l_space_before;
0081
0082 for (int j = 0; j < repeat; j++) {
0083 string l_name = Form("layer%d", l_num);
0084 double l_thickness = layering.layer(l_num - 1)->thickness();
0085 double l_dim_x = tan_hphi* l_pos_z;
0086 l_pos_z += l_thickness;
0087
0088 Position l_pos(0, 0, l_pos_z - l_thickness/2.);
0089 double l_trd_x1 = l_dim_x;
0090 double l_trd_x2 = l_dim_x + l_thickness*tan_hphi;
0091 double l_trd_y1 = l_dim_y;
0092 double l_trd_y2 = l_trd_y1;
0093 double l_trd_z = l_thickness/2;
0094 Trapezoid l_shape(l_trd_x1, l_trd_x2, l_trd_y1, l_trd_y2, l_trd_z);
0095 Volume l_vol(l_name, l_shape, air);
0096 DetElement layer(stave_det, l_name, det_id);
0097
0098
0099 int s_num = 1;
0100 double s_pos_z = -(l_thickness / 2.);
0101 for(xml_coll_t si(x_layer,_U(slice)); si; ++si) {
0102 xml_comp_t x_slice = si;
0103 string s_name = Form("slice%d", s_num);
0104 double s_thick = x_slice.thickness();
0105 double s_trd_x1 = l_dim_x + (s_pos_z + l_thickness/2)*tan_hphi;
0106 double s_trd_x2 = l_dim_x + (s_pos_z + l_thickness/2 + s_thick)*tan_hphi;
0107 double s_trd_y1 = l_trd_y1;
0108 double s_trd_y2 = s_trd_y1;
0109 double s_trd_z = s_thick/2.;
0110 Trapezoid s_shape(s_trd_x1, s_trd_x2, s_trd_y1, s_trd_y2, s_trd_z);
0111 Volume s_vol(s_name, s_shape, desc.material(x_slice.materialStr()));
0112 DetElement slice(layer, s_name, det_id);
0113
0114
0115 if (x_slice.hasChild(_Unicode(fiber))) {
0116 buildFibers(desc, sens, s_vol, x_slice.child(_Unicode(fiber)), {s_trd_x1, s_thick, l_dim_y, hphi});
0117 }
0118
0119 if ( x_slice.isSensitive() ) {
0120 s_vol.setSensitiveDetector(sens);
0121 }
0122 s_vol.setAttributes(desc, x_slice.regionStr(), x_slice.limitsStr(), x_slice.visStr());
0123
0124
0125 PlacedVolume slice_phv = l_vol.placeVolume(s_vol, Position(0, 0, s_pos_z + s_thick/2));
0126 slice_phv.addPhysVolID("slice", s_num);
0127 slice.setPlacement(slice_phv);
0128
0129 s_pos_z += s_thick;
0130 ++s_num;
0131 }
0132
0133
0134 l_vol.setAttributes(desc, x_layer.regionStr(), x_layer.limitsStr(), x_layer.visStr());
0135
0136 PlacedVolume layer_phv = mod_vol.placeVolume(l_vol, l_pos);
0137 layer_phv.addPhysVolID("layer", l_num);
0138 layer.setPlacement(layer_phv);
0139
0140 if (j < repeat - 1) {
0141 l_pos_z += l_space_between;
0142 }
0143 ++l_num;
0144 }
0145 }
0146 }
0147
0148 double phi = M_PI / nsides;
0149
0150 for (int i = 0; i < nsides; i++, phi -= dphi) {
0151
0152 Transform3D tr(RotationZYX(0, phi, M_PI*0.5), Translation3D(0, 0, 0));
0153 PlacedVolume pv = envelope.placeVolume(mod_vol, tr);
0154 pv.addPhysVolID("module", i + 1);
0155 DetElement sd = (i == 0) ? stave_det : stave_det.clone(Form("stave%d", i));
0156 sd.setPlacement(pv);
0157 sdet.add(sd);
0158 }
0159
0160
0161 if (x_det.hasChild(_U(staves))) {
0162 xml_comp_t x_staves = x_det.staves();
0163 mod_vol.setVisAttributes(desc.visAttributes(x_staves.visStr()));
0164 if (x_staves.hasChild(_U(support))) {
0165 buildSupport(desc, mod_vol, x_staves.child(_U(support)), {inner_r, l_pos_z, x_dim.z(), hphi});
0166 }
0167 }
0168
0169
0170 envelope.setAttributes(desc, x_det.regionStr(), x_det.limitsStr(), x_det.visStr());
0171 return sdet;
0172 }
0173
0174 void buildFibers(Detector& desc, SensitiveDetector &sens, Volume &s_vol, xml_comp_t x_fiber,
0175 const std::tuple<double, double, double, double> &dimensions)
0176 {
0177 auto [s_trd_x1, s_thick, s_length, hphi] = dimensions;
0178 double f_radius = getAttrOrDefault(x_fiber, _U(radius), 0.1 * cm);
0179 double f_spacing_x = getAttrOrDefault(x_fiber, _Unicode(spacing_x), 0.122 * cm);
0180 double f_spacing_z = getAttrOrDefault(x_fiber, _Unicode(spacing_z), 0.134 * cm);
0181 std::string f_id_grid = getAttrOrDefault(x_fiber, _Unicode(identifier_grid), "grid");
0182 std::string f_id_fiber = getAttrOrDefault(x_fiber, _Unicode(identifier_fiber), "fiber");
0183
0184
0185
0186
0187
0188
0189
0190
0191
0192 auto grid_div = getNdivisions(s_trd_x1, s_thick, 2.0*cm, 2.0*cm);
0193
0194 auto grid_vtx = gridPoints(grid_div.first, grid_div.second, s_trd_x1, s_thick, hphi);
0195 Tube f_tube(0, f_radius, s_length);
0196 Volume f_vol("fiber_vol", f_tube, desc.material(x_fiber.materialStr()));
0197
0198 vector<int> f_id_count(grid_div.first*grid_div.second, 0);
0199 auto f_pos = fiberPositions(f_radius, f_spacing_x, f_spacing_z, s_trd_x1, s_thick, hphi);
0200
0201 for (size_t il = 0; il < f_pos.size(); ++il) {
0202 auto &line = f_pos[il];
0203 if (line.empty()) {
0204 continue;
0205 }
0206 double l_pos_y = line.front().y();
0207
0208 Assembly lfibers(Form("fiber_array_line_%lu", il));
0209 for (auto &p : line) {
0210 int f_grid_id = -1;
0211 int f_id = -1;
0212
0213 for (auto &poly_vtx : grid_vtx) {
0214 if (p.y() != l_pos_y) {
0215 std::cerr << Form("Expected the same y position from a same line: %.2f, but got %.2f", l_pos_y, p.y())
0216 << std::endl;
0217 continue;
0218 }
0219 auto [grid_id, vtx_a, vtx_b, vtx_c, vtx_d] = poly_vtx;
0220 double poly_x[4] = {vtx_a.x(), vtx_b.x(), vtx_c.x(), vtx_d.x()};
0221 double poly_y[4] = {vtx_a.y(), vtx_b.y(), vtx_c.y(), vtx_d.y()};
0222 double f_xy[2] = {p.x(), p.y()};
0223
0224 TGeoPolygon poly(4);
0225 poly.SetXY(poly_x, poly_y);
0226 poly.FinishPolygon();
0227
0228 if(poly.Contains(f_xy)) {
0229 f_grid_id = grid_id;
0230 f_id = f_id_count[grid_id];
0231 f_id_count[grid_id]++;
0232 }
0233 }
0234
0235 if ( x_fiber.isSensitive() ) {
0236 f_vol.setSensitiveDetector(sens);
0237 }
0238 f_vol.setAttributes(desc, x_fiber.regionStr(), x_fiber.limitsStr(), x_fiber.visStr());
0239
0240
0241
0242
0243 PlacedVolume fiber_phv = lfibers.placeVolume(f_vol, Position(p.x(), 0., 0.));
0244 fiber_phv.addPhysVolID(f_id_grid, f_grid_id + 1).addPhysVolID(f_id_fiber, f_id + 1);
0245 }
0246 lfibers.ptr()->Voxelize("");
0247 Transform3D l_tr(RotationZYX(0,0,M_PI*0.5),Position(0., 0, l_pos_y));
0248 s_vol.placeVolume(lfibers, l_tr);
0249 }
0250 }
0251
0252
0253
0254 void buildSupport(Detector& desc, Volume &mod_vol, xml_comp_t x_support,
0255 const std::tuple<double, double, double, double> &dimensions)
0256 {
0257 auto [inner_r, l_pos_z, stave_length, hphi] = dimensions;
0258
0259 double support_thickness = getAttrOrDefault(x_support, _Unicode(thickness), 5. * cm);
0260 double beam_thickness = getAttrOrDefault(x_support, _Unicode(beam_thickness), support_thickness/4.);
0261
0262 if (beam_thickness > support_thickness/3.) {
0263 std::cerr << Form("beam_thickness (%.2f) cannot be greater than support_thickness/3 (%.2f), shrink it to fit",
0264 beam_thickness, support_thickness/3.) << std::endl;
0265 beam_thickness = support_thickness/3.;
0266 }
0267 double trd_x1_support = std::tan(hphi) * l_pos_z;
0268 double trd_x2_support = std::tan(hphi) * (l_pos_z + support_thickness);
0269 double trd_y = stave_length / 2.;
0270 Assembly env_vol ("support_envelope");
0271
0272 double grid_size = getAttrOrDefault(x_support, _Unicode(grid_size), 25. * cm);
0273 int n_cross_supports = std::floor(trd_y - beam_thickness)/grid_size;
0274
0275
0276 int n_beams = getAttrOrDefault(x_support, _Unicode(n_beams), 3);;
0277 double beam_width = 2. * trd_x1_support / (n_beams + 1);
0278 double beam_gap = getAttrOrDefault(x_support, _Unicode(beam_gap), 3.*cm);
0279
0280
0281 double beam_space_x = beam_width + beam_gap;
0282 double beam_space_z = support_thickness - beam_thickness;
0283 double cross_thickness = support_thickness - beam_thickness;
0284 double beam_pos_z = beam_thickness / 2.;
0285 double beam_center_z = support_thickness / 2. - beam_pos_z;
0286
0287 Box beam_vert_s(beam_thickness / 2., trd_y, cross_thickness / 2.);
0288 Box beam_hori_s(beam_width / 2., trd_y, beam_thickness / 2.);
0289 UnionSolid T_beam_s(beam_hori_s, beam_vert_s, Position(0., 0., support_thickness / 2.));
0290 Volume H_beam_vol("H_beam", T_beam_s, desc.material(x_support.materialStr()));
0291 H_beam_vol.setVisAttributes(desc, x_support.visStr());
0292
0293 double beam_start_x = - (n_beams - 1) * (beam_width + beam_gap) / 2.;
0294 for (int i = 0; i < n_beams; ++i) {
0295 Position beam_pos(beam_start_x + i * (beam_width + beam_gap), 0., - support_thickness / 2. + beam_pos_z);
0296 env_vol.placeVolume(H_beam_vol, beam_pos);
0297 }
0298
0299
0300 double cross_x = beam_space_x - beam_thickness;
0301 Box cross_s(cross_x / 2., beam_thickness / 2., cross_thickness / 2.);
0302 Volume cross_vol("cross_center_beam", cross_s, desc.material(x_support.materialStr()));
0303 cross_vol.setVisAttributes(desc, x_support.visStr());
0304 for (int i = 0; i < n_beams - 1; ++i) {
0305 env_vol.placeVolume(cross_vol, Position(beam_start_x + beam_space_x * (i + 0.5), 0., beam_pos_z));
0306 for (int j = 1; j < n_cross_supports; j++) {
0307 env_vol.placeVolume(cross_vol, Position(beam_start_x + beam_space_x * (i + 0.5), -j * grid_size, beam_pos_z));
0308 env_vol.placeVolume(cross_vol, Position(beam_start_x + beam_space_x * (i + 0.5), j * grid_size, beam_pos_z));
0309 }
0310 }
0311
0312
0313
0314 double cross_edge_x = trd_x1_support + beam_start_x - beam_thickness / 2.;
0315 double cross_trd_x1 = cross_edge_x + std::tan(hphi) * beam_thickness;
0316 double cross_trd_x2 = cross_trd_x1 + 2.* std::tan(hphi) * cross_thickness;
0317 double edge_pos_x = beam_start_x - cross_trd_x1 / 2. - beam_thickness / 2;
0318 Trapezoid cross_s2_trd (cross_trd_x1 / 2., cross_trd_x2 / 2.,
0319 beam_thickness / 2., beam_thickness / 2., cross_thickness / 2.);
0320 Box cross_s2_box ((cross_trd_x2 - cross_trd_x1)/4., beam_thickness / 2., cross_thickness / 2.);
0321 SubtractionSolid cross_s2(cross_s2_trd, cross_s2_box, Position((cross_trd_x2 + cross_trd_x1)/4., 0., 0.));
0322 Volume cross_vol2("cross_edge_beam", cross_s2, desc.material(x_support.materialStr()));
0323 cross_vol2.setVisAttributes(desc, x_support.visStr());
0324 env_vol.placeVolume(cross_vol2, Position(edge_pos_x, 0., beam_pos_z));
0325 env_vol.placeVolume(cross_vol2, Transform3D(Translation3D(-edge_pos_x, 0., beam_pos_z) * RotationZ(M_PI)));
0326 for (int j = 1; j < n_cross_supports; j++) {
0327 env_vol.placeVolume(cross_vol2, Position(edge_pos_x, -j * grid_size, beam_pos_z));
0328 env_vol.placeVolume(cross_vol2, Position(edge_pos_x, j * grid_size, beam_pos_z));
0329 env_vol.placeVolume(cross_vol2, Transform3D(Translation3D(-edge_pos_x, -j * grid_size, beam_pos_z) * RotationZ(M_PI)));
0330 env_vol.placeVolume(cross_vol2, Transform3D(Translation3D(-edge_pos_x, j * grid_size, beam_pos_z) * RotationZ(M_PI)));
0331 }
0332
0333 mod_vol.placeVolume(env_vol, Position(0.0, 0.0, l_pos_z + support_thickness/2.));
0334 }
0335
0336
0337 vector<vector<Point>> fiberPositions(double radius, double x_spacing, double z_spacing, double x, double z, double phi,
0338 double spacing_tol)
0339 {
0340
0341
0342
0343
0344
0345
0346 vector<vector<Point>> positions;
0347 int z_layers = floor((z / 2 - radius - spacing_tol) / z_spacing);
0348
0349 double z_pos = 0.;
0350 double x_pos = 0.;
0351
0352 for (int l = -z_layers; l < z_layers + 1; l++) {
0353 vector<Point> xline;
0354 z_pos = l * z_spacing;
0355 double x_max = x + (z / 2. + z_pos) * tan(phi) - spacing_tol;
0356 (l % 2 == 0) ? x_pos = 0. : x_pos = x_spacing / 2;
0357
0358 while (x_pos < (x_max - radius)) {
0359 xline.push_back(Point(x_pos, z_pos));
0360 if (x_pos != 0.)
0361 xline.push_back(Point(-x_pos, z_pos));
0362 x_pos += x_spacing;
0363 }
0364
0365 sort(xline.begin(), xline.end(), [](const Point& p1, const Point& p2) { return p1.x() < p2.x(); });
0366 positions.emplace_back(std::move(xline));
0367 }
0368 return positions;
0369 }
0370
0371
0372 std::pair<int, int> getNdivisions(double x, double z, double dx, double dz)
0373 {
0374
0375
0376
0377
0378 double SiPMsize = 13.0 * mm;
0379 double grid_min = SiPMsize + 3.0 * mm;
0380
0381 if (dz < grid_min) {
0382 dz = grid_min;
0383 }
0384
0385 if (dx < grid_min) {
0386 dx = grid_min;
0387 }
0388
0389 int nfit_cells_z = floor(z / dz);
0390 int n_cells_z = nfit_cells_z;
0391
0392 if (nfit_cells_z == 0)
0393 n_cells_z++;
0394
0395 int nfit_cells_x = floor((2 * x) / dx);
0396 int n_cells_x = nfit_cells_x;
0397
0398 if (nfit_cells_x == 0)
0399 n_cells_x++;
0400
0401 return std::make_pair(n_cells_x, n_cells_z);
0402 }
0403
0404
0405 vector<tuple<int, Point, Point, Point, Point>> gridPoints(int div_x, int div_z, double x, double z, double phi)
0406 {
0407
0408
0409 double dz = z / div_z;
0410
0411 std::vector<std::tuple<int, Point, Point, Point, Point>> points;
0412
0413 for (int iz = 0; iz < div_z + 1; iz++) {
0414 for (int ix = 0; ix < div_x + 1; ix++) {
0415 double A_z = -z / 2 + iz * dz;
0416 double B_z = -z / 2 + (iz + 1) * dz;
0417
0418 double len_x_for_z = 2 * (x + iz * dz * tan(phi));
0419 double len_x_for_z_plus_1 = 2 * (x + (iz + 1) * dz * tan(phi));
0420
0421 double dx_for_z = len_x_for_z / div_x;
0422 double dx_for_z_plus_1 = len_x_for_z_plus_1 / div_x;
0423
0424 double A_x = -len_x_for_z / 2. + ix * dx_for_z;
0425 double B_x = -len_x_for_z_plus_1 / 2. + ix * dx_for_z_plus_1;
0426
0427 double C_z = B_z;
0428 double D_z = A_z;
0429 double C_x = B_x + dx_for_z_plus_1;
0430 double D_x = A_x + dx_for_z;
0431
0432 int id = ix + div_x * iz;
0433
0434 auto A = Point(A_x, A_z);
0435 auto B = Point(B_x, B_z);
0436 auto C = Point(C_x, C_z);
0437 auto D = Point(D_x, D_z);
0438
0439
0440 points.push_back(make_tuple(id, A, B, C, D));
0441 }
0442 }
0443
0444 return points;
0445 }
0446
0447 DECLARE_DETELEMENT(athena_EcalBarrelInterlayers, create_detector)
0448