epic/src/DIRC_geo.cpp

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 - 2026 Wouter Deconinck, Dmitry Romanov, Bill Llope
0003
0004 #include "DD4hep/DetFactoryHelper.h"
0005 #include "DD4hep/OpticalSurfaces.h"
0006 #include "DD4hep/Printout.h"
0007 #include "DDRec/DetectorData.h"
0008 #include "DDRec/Surface.h"
0009 #include <XML/Helper.h>
0010
0011 //////////////////////////////////
0012 // Central Barrel DIRC
0013 //////////////////////////////////
0014
0015 using namespace std;
0016 using namespace dd4hep;
0017
0018 static dd4hep::Trap MakeTrap(const std::string& pName, double pZ, double pY, double pX,
0019                              double pLTX);
0020
0021 //--------------------------------------------------------------------------------
0022 //
0023 static Ref_t createDetector(Detector& desc, xml_h e, SensitiveDetector sens) {
0024   xml_det_t xml_det = e;
0025
0026 //#define VERBOSE_DIRCgeo
0027 #ifdef VERBOSE_DIRCgeo
0028   // TEMPORARILY INCREASE VERBOSITY level for debugging purposes
0029   PrintLevel priorPrintLevel = printLevel();
0030   setPrintLevel(DEBUG);
0031 #endif
0032
0033   //---- steering section
0034   //
0035   double CADD_barplaneR_exp       = 765.875 * mm;  // CADD, Kris Cleveland, 11/14/2025
0036   double CADD_mirrorthickness_exp = 0.3 * cm;      // CADD, Kris Cleveland, 11/14/2025
0037   double FAKE_dirc_pos_z_exp      = -44.7525 * cm; //
0038   double CADD_gluethickness_exp   = 0.025 * mm;    // CADD, Kris Cleveland, 11/14/2025
0039   double CADD_barheight_exp       = 17.25 * mm;    // CADD, Kris Cleveland, 11/14/2025
0040   double CADD_barwidth_exp        = 35.00 * mm;    //
0041   double FAKE_barlength_exp =
0042       1140.000 * mm; // fake bar length to make 4 bars have correct CADD length...
0043   //double CADD_motherZshift            =     3.075 * mm - 0.66*cm;     // to put mirror Zmax at CADD_mirrorZmax
0044   //double FAKE_barlength                       =  1136.7   * mm;               // ORIG VALUE
0045   double barsZmin_expected   = -2729.075 * mm; //
0046   double barsZmax_expected   = 1831.025 * mm;  //
0047   double barsZlen_expected   = barsZmax_expected - barsZmin_expected;
0048   double mirrorZmax_expected = 1834.025 * mm; // CADD, Kris Cleveland, 11/14/2025
0049   double RailRadius          = 77.05 * cm;    // ORIG VALUE, needed for rails to avoid conflicts...
0050   printout(DEBUG, "DIRC_geo", "---- hpDIRC geometry. All numbers below in cm ----");
0051   //
0052   //---- end steering section
0053
0054   //---- Detector element
0055   string det_name = xml_det.nameStr();
0056   int det_id      = xml_det.id();
0057   DetElement det(det_name, det_id);
0058   printout(DEBUG, "DIRC_geo", "det_name        = %s", det_name.c_str());
0059   printout(DEBUG, "DIRC_geo", "det_id          = %d", det_id);
0060
0061   //---- Detector dimension, position, rotation
0062   xml_dim_t dirc_dim = xml_det.dimensions();
0063   xml_dim_t dirc_pos = xml_det.position();
0064   double det_rmin    = dirc_dim.rmin();
0065   double det_rmax    = dirc_dim.rmax();
0066   double det_ravg    = (det_rmin + det_rmax) / 2.;
0067   printout(DEBUG, "DIRC_geo", "det_rmin           = %12.4f", det_rmin);
0068   printout(DEBUG, "DIRC_geo", "det_rmax           = %12.4f", det_rmax);
0069   printout(DEBUG, "DIRC_geo", "det_ravg           = %12.4f", det_ravg);
0070   printout(DEBUG, "DIRC_geo", "CADD_barplaneR_exp = %12.4f", CADD_barplaneR_exp);
0071   printout(DEBUG, "DIRC_geo", "RailRadius         = %12.4f", RailRadius);
0072   //
0073   double dirc_pos_z = dirc_pos.z();
0074   printout(DEBUG, "DIRC_geo", "dirc_pos_z      = %12.4f", dirc_pos_z);
0075   printout(DEBUG, "DIRC_geo", "dirc_pos_z exp  = %12.4f", FAKE_dirc_pos_z_exp);
0076
0077   //---- Detector type
0078   sens.setType("tracker");
0079
0080   //---- Entire DIRC assembly
0081   Assembly det_volume("DIRC");
0082   det_volume.setVisAttributes(desc.visAttributes(xml_det.visStr()));
0083   Transform3D det_tr(RotationY(0), Position(0.0, 0.0, dirc_pos_z));
0084   det.setPlacement(
0085       desc.pickMotherVolume(det).placeVolume(det_volume, det_tr).addPhysVolID("system", det_id));
0086
0087   //---- Assembly dirc_module holds all the hpDIRC
0088   //
0089   xml_comp_t xml_module = xml_det.child(_U(module));
0090   Assembly dirc_module("DIRCModule");
0091   dirc_module.setVisAttributes(desc.visAttributes(xml_module.visStr()));
0092
0093   //---- define bar_vol
0094   xml_comp_t xml_bar = xml_module.child(_Unicode(bar));
0095   double bar_height  = xml_bar.height();
0096   double bar_width   = xml_bar.width();
0097   double bar_length  = xml_bar.length();
0098   printout(DEBUG, "DIRC_geo", "bar_height         = %12.4f", bar_height);
0099   printout(DEBUG, "DIRC_geo", "CADD_barheight_exp = %12.4f", CADD_barheight_exp);
0100   printout(DEBUG, "DIRC_geo", "bar_width          = %12.4f", bar_width);
0101   printout(DEBUG, "DIRC_geo", "CADD_barwidth_exp  = %12.4f", CADD_barwidth_exp);
0102   printout(DEBUG, "DIRC_geo", "bar_length         = %12.4f", bar_length);
0103   printout(DEBUG, "DIRC_geo", "FAKE_barlength_exp = %12.4f", FAKE_barlength_exp);
0104   //
0105   Box bar_box("bar_box", bar_height / 2, bar_width / 2,
0106               bar_length / 2); // the box volume for a bar, not a "bar box"
0107   Volume bar_vol("bar_vol", bar_box, desc.material(xml_bar.materialStr()));
0108   bar_vol.setVisAttributes(desc.visAttributes(xml_bar.visStr()));
0109
0110   //---- define glue_vol
0111   xml_comp_t xml_glue   = xml_module.child(_Unicode(glue));
0112   double glue_thickness = xml_glue.thickness();
0113   Box glue_box("glue_box", bar_height / 2, bar_width / 2, glue_thickness / 2);
0114   Volume glue_vol("glue_vol", glue_box, desc.material(xml_glue.materialStr()));
0115   glue_vol.setVisAttributes(desc.visAttributes(xml_glue.visStr()));
0116   printout(DEBUG, "DIRC_geo", "glue_thickness         = %12.4f", glue_thickness);
0117   printout(DEBUG, "DIRC_geo", "CADD_gluethickness_exp = %12.4f", CADD_gluethickness_exp);
0118
0119   //---- barbox contents variables
0120   auto bar_repeat_y    = xml_bar.attr<int>(_Unicode(repeat_y));
0121   auto bar_repeat_z    = xml_bar.attr<int>(_Unicode(repeat_z));
0122   auto bar_gap         = xml_bar.gap();
0123   auto bar_assm_width  = bar_width * bar_repeat_y + bar_gap * (bar_repeat_y - 1);
0124   auto bar_assm_length = (bar_length + glue_thickness) * bar_repeat_z;
0125   printout(DEBUG, "DIRC_geo", "bar_repeat_y        = %8d ", bar_repeat_y);
0126   printout(DEBUG, "DIRC_geo", "bar_repeat_z        = %8d ", bar_repeat_z);
0127   printout(DEBUG, "DIRC_geo", "bar_gap             = %12.4f ", bar_gap);
0128   printout(DEBUG, "DIRC_geo", "bar_assm_width      = %12.4f ", bar_assm_width);
0129   printout(DEBUG, "DIRC_geo", "bar_assm_length     = %12.4f ", bar_assm_length);
0130
0131   //---- Mirror construction
0132   xml_comp_t xml_mirror = xml_module.child(_Unicode(mirror));
0133   auto mirror_width     = xml_mirror.width();
0134   auto mirror_height    = xml_mirror.height();
0135   auto mirror_thickness = xml_mirror.thickness();
0136   printout(DEBUG, "DIRC_geo", "mirror_width             = %12.4f", mirror_width);
0137   printout(DEBUG, "DIRC_geo", "mirror_height            = %12.4f", mirror_height);
0138   printout(DEBUG, "DIRC_geo", "mirror_thickness         = %12.4f", mirror_thickness);
0139   printout(DEBUG, "DIRC_geo", "CADD_mirrorthickness_exp = %12.4f", CADD_mirrorthickness_exp);
0140   //
0141   Box mirror_box("mirror_box", mirror_height / 2, mirror_width / 2, mirror_thickness / 2);
0142   Volume mirror_vol("mirror_vol", mirror_box, desc.material(xml_mirror.materialStr()));
0143   mirror_vol.setVisAttributes(desc.visAttributes(xml_mirror.visStr()));
0144   //
0145   //---- Mirror optical surface
0146   auto surfMgr = desc.surfaceManager();
0147   auto surf    = surfMgr.opticalSurface("DIRC_MirrorOpticalSurface");
0148   SkinSurface skin(desc, det, Form("dirc_mirror_optical_surface"), surf, mirror_vol);
0149   skin.isValid();
0150
0151   //---- Define Envelope_box_vol that holds 4*10bars+glue+mirror
0152   //---- Place Envelope_box_vol inside dirc_module
0153   //    ...note Envelope_box_vol is long enough to contain the bars and mirror
0154   //    ...dirc_module is placed inside Envelope_box_vol with z-offset that is half mirror_thickness.
0155   //
0156   auto envbox_xsize = (mirror_height + 1 * mm) / 2;
0157   auto envbox_ysize = (mirror_width + 1 * mm) / 2;
0158   auto envbox_zsize = 0.5 * bar_repeat_z * (bar_length + glue_thickness) + 0.5 * mirror_thickness;
0159   Box Envelope_box("Envelope_box", envbox_xsize, envbox_ysize, envbox_zsize);
0160   Volume Envelope_box_vol("Envelope_box_vol", Envelope_box, desc.material("AirOptical"));
0161   dirc_module.placeVolume(Envelope_box_vol, Position(0, 0, 0.5 * mirror_thickness));
0162   printout(DEBUG, "DIRC_geo", "envbox_xsize    = %12.4f ", envbox_xsize);
0163   printout(DEBUG, "DIRC_geo", "envbox_ysize    = %12.4f ", envbox_ysize);
0164   printout(DEBUG, "DIRC_geo", "envbox_zsize    = %12.4f ", envbox_zsize);
0165   printout(DEBUG, "DIRC_geo", "envbox_xpos     = %12.4f ", 0);
0166   printout(DEBUG, "DIRC_geo", "envbox_ypos     = %12.4f ", 0);
0167   printout(DEBUG, "DIRC_geo", "envbox_zpos LAB = %12.4f ", 0.5 * mirror_thickness + dirc_pos_z);
0168
0169   //---- place glue_vol and bar_vol into Envelope_box_vol
0170   //
0171   //            Note z_index = 0 is the high-Z bar! (pos eta)
0172   //            Note z_index = 3 is the  low-Z bar! (neg eta)
0173   //
0174   //            Note: bar's x (depth) position is same as that of dirc mother boxes.
0175   //
0176   double barZmax = 0;
0177   double barZmin = 0;
0178   for (int y_index = 0; y_index < bar_repeat_y; y_index++) {
0179     double y = 0.5 * bar_assm_width - 0.5 * bar_width - (bar_width + bar_gap) * y_index;
0180     for (int z_index = 0; z_index < bar_repeat_z; z_index++) {
0181       double z = 0.5 * bar_assm_length - 0.5 * mirror_thickness - 0.5 * bar_length -
0182                  (bar_length + glue_thickness) * z_index;
0183       Envelope_box_vol.placeVolume(glue_vol,
0184                                    Position(0, y, z - 0.5 * (bar_length + glue_thickness)));
0185       Envelope_box_vol.placeVolume(bar_vol, Position(0, y, z))
0186           .addPhysVolID("section", z_index)
0187           .addPhysVolID("bar", y_index);
0188       //
0189       if (y_index == 4 && z_index == 0) {
0190         barZmax = z + bar_length / 2.;
0191       }
0192       if (y_index == 4 && z_index == 3) {
0193         barZmin = z - bar_length / 2. - glue_thickness;
0194       }
0195       printout(DEBUG, "DIRC_geo",
0196                "iy,iz: %3d,%3d \t bar y,z: %10.4f,%10.4f \t glue y,z: %10.4f,%10.4f", y_index,
0197                z_index, y, z + dirc_pos_z, y, z + dirc_pos_z - 0.5 * (bar_length + glue_thickness));
0198       //
0199     }
0200   }
0201   printout(DEBUG, "DIRC_geo", "barZmax LAB          = %12.4f ", barZmax + dirc_pos_z);
0202   printout(DEBUG, "DIRC_geo", "barZmax LAB expected = %12.4f ", barsZmax_expected);
0203   printout(DEBUG, "DIRC_geo", "barZmin LAB          = %12.4f ", barZmin + dirc_pos_z);
0204   printout(DEBUG, "DIRC_geo", "barZmin LAB expected = %12.4f ", barsZmin_expected);
0205   printout(DEBUG, "DIRC_geo", "barZlen LAB          = %12.4f ", barZmax - barZmin);
0206   printout(DEBUG, "DIRC_geo", "barZlen LAB expected = %12.4f ", barsZlen_expected);
0207
0208   //---- Place mirror
0209   Envelope_box_vol.placeVolume(mirror_vol, Position(0, 0, 0.5 * bar_assm_length));
0210   printout(DEBUG, "DIRC_geo", "mirror_zpos LAB          = %12.4f ",
0211            0.5 * bar_assm_length + dirc_pos_z);
0212   printout(DEBUG, "DIRC_geo", "mirror Zmax LAB          = %12.4f ",
0213            0.5 * bar_assm_length + dirc_pos_z + mirror_thickness / 2.);
0214   printout(DEBUG, "DIRC_geo", "mirror Zmax LAB expected = %12.4f ", mirrorZmax_expected);
0215
0216   //---- Prism variables
0217   xml_comp_t xml_prism    = xml_module.child(_Unicode(prism));
0218   double prism_angle      = xml_prism.angle();
0219   double prism_width      = xml_prism.width();
0220   double prism_length     = xml_prism.length();
0221   double prism_short_edge = getAttrOrDefault(xml_prism, _Unicode(short_edge), 50 * mm);
0222   double prism_long_edge  = prism_short_edge + prism_length * tan(prism_angle);
0223   printout(DEBUG, "DIRC_geo", "prism_angle       = %12.4f ", prism_angle);
0224   printout(DEBUG, "DIRC_geo", "prism_width       = %12.4f ", prism_width);
0225   printout(DEBUG, "DIRC_geo", "prism_length      = %12.4f ", prism_length);
0226   printout(DEBUG, "DIRC_geo", "prism_short_edge  = %12.4f ", prism_short_edge);
0227   printout(DEBUG, "DIRC_geo", "prism_long_edge   = %12.4f ", prism_long_edge);
0228
0229   //---- Lens variables
0230   xml_comp_t xml_lens   = xml_module.child(_Unicode(lens));
0231   double lens_shift     = getAttrOrDefault(xml_lens, _Unicode(shift), 0 * mm);
0232   double lens_width     = getAttrOrDefault(xml_lens, _Unicode(width), 35 * mm);
0233   double lens_thickness = getAttrOrDefault(xml_lens, _Unicode(thickness), 12 * mm);
0234   double lens_r1        = getAttrOrDefault(xml_lens, _Unicode(r1), 62 * mm);
0235   double lens_r2        = getAttrOrDefault(xml_lens, _Unicode(r2), 36 * mm);
0236
0237   //---- Lens construction: 3-layer spherical lens ---
0238   double lens_radius        = sqrt(lens_width * lens_width / 4. + bar_height * bar_height / 4.);
0239   double lens_min_thickness = 2.0 * mm;
0240   double ztrans1 = -lens_thickness / 2. - sqrt(lens_r1 * lens_r1 - lens_radius * lens_radius) +
0241                    lens_min_thickness;
0242   double ztrans2 = -lens_thickness / 2. - sqrt(lens_r2 * lens_r2 - lens_radius * lens_radius) +
0243                    lens_min_thickness * 2;
0244   Box lens_symm_box("lens_symm_box", 0.5 * prism_short_edge, 0.5 * lens_width,
0245                     0.5 * lens_thickness);
0246   Volume Envelope_lens_vol("Envelope_lens_vol", lens_symm_box, desc.material("AirOptical"));
0247   Tube lens_symm_tube(0, lens_radius, 0.5 * lens_thickness);
0248   Sphere lens_sphere1(0, lens_r1);
0249   Sphere lens_sphere2(0, lens_r2);
0250   IntersectionSolid lens_box("lens_box", lens_symm_box, lens_symm_box,
0251                              Position(0, 0, -lens_min_thickness * 2));
0252   IntersectionSolid lens_tube("lens_tube", lens_symm_tube, lens_symm_box,
0253                               Position(0, 0, lens_min_thickness * 2));
0254   UnionSolid lens_box_tube("lens_box_tube", lens_box, lens_tube);
0255   IntersectionSolid lens_layer1_solid("lens_layer1_solid", lens_box_tube, lens_sphere1,
0256                                       Position(0, 0, -ztrans1));
0257   SubtractionSolid lens_layer23_solid("lens_layer23_solid", lens_box_tube, lens_sphere1,
0258                                       Position(0, 0, -ztrans1));
0259   IntersectionSolid lens_layer2_solid("lens_layer2_solid", lens_layer23_solid, lens_sphere2,
0260                                       Position(0, 0, -ztrans2));
0261   SubtractionSolid lens_layer3_solid("lens_layer3_solid", lens_layer23_solid, lens_sphere2,
0262                                      Position(0, 0, -ztrans2));
0263   Volume lens_layer1_vol("lens_layer1_vol", lens_layer1_solid,
0264                          desc.material(xml_lens.attr<std::string>(_Unicode(material1))));
0265   Volume lens_layer2_vol("lens_layer2_vol", lens_layer2_solid,
0266                          desc.material(xml_lens.attr<std::string>(_Unicode(material2))));
0267   Volume lens_layer3_vol("lens_layer3_vol", lens_layer3_solid,
0268                          desc.material(xml_lens.attr<std::string>(_Unicode(material3))));
0269   lens_layer1_vol.setVisAttributes(desc.visAttributes(xml_lens.attr<std::string>(_Unicode(vis1))));
0270   lens_layer2_vol.setVisAttributes(desc.visAttributes(xml_lens.attr<std::string>(_Unicode(vis2))));
0271   lens_layer3_vol.setVisAttributes(desc.visAttributes(xml_lens.attr<std::string>(_Unicode(vis3))));
0272   //
0273   double lens_position_x = lens_shift;
0274   double lens_position_z = -0.5 * (bar_assm_length + lens_thickness);
0275   //std::cout<<"lens_position_x     = "<<lens_position_x <<std::endl;
0276   //std::cout<<"lens_position_z     = "<<lens_position_z <<std::endl;
0277   printout(DEBUG, "DIRC_geo", "lens_position_x LAB = %12.4f ", lens_position_x + det_ravg);
0278   printout(DEBUG, "DIRC_geo", "lens_position_z LAB = %12.4f ", lens_position_z + dirc_pos_z);
0279   //
0280   for (int y_index = 0; y_index < bar_repeat_y; y_index++) {
0281     double lens_position_y = y_index * lens_width - 0.5 * (prism_width - lens_width);
0282     Position lens_position(lens_position_x, lens_position_y, lens_position_z);
0283     dirc_module.placeVolume(Envelope_lens_vol, lens_position);
0284   }
0285   Envelope_lens_vol.placeVolume(lens_layer1_vol);
0286   Envelope_lens_vol.placeVolume(lens_layer2_vol);
0287   Envelope_lens_vol.placeVolume(lens_layer3_vol);
0288
0289   //---- Prism construction
0290   Trap prism_trap =
0291       MakeTrap("prism_trap", prism_width, prism_length, prism_long_edge, prism_short_edge);
0292   Volume prism_vol("prism_vol", prism_trap, desc.material(xml_prism.materialStr()));
0293   prism_vol.setVisAttributes(desc.visAttributes(xml_prism.visStr()));
0294   //
0295   double prism_position_x =
0296       (prism_long_edge + prism_short_edge) / 4. - 0.5 * prism_short_edge + lens_shift;
0297   double prism_position_z = -0.5 * (bar_assm_length + prism_length) - lens_thickness;
0298   RotationX prism_rotation(M_PI / 2.);
0299   Position prism_position(prism_position_x, 0, prism_position_z);
0300   printout(DEBUG, "DIRC_geo", "prism_position_x LAB = %12.4f ", prism_position_x + det_ravg);
0301   printout(DEBUG, "DIRC_geo", "prism_position_z LAB = %12.4f ", prism_position_z + dirc_pos_z);
0302
0303   //---- Envelope for prism + mcp
0304   //    Envelope_trap_vol contains prism_vol and mcp_vol
0305   //
0306   double Envelope_trap_width      = prism_width + 1 * mm;
0307   double Envelope_trap_length     = prism_length + 1 * mm; // mcp thickness is 1 mm
0308   double Envelope_trap_short_edge = prism_short_edge + 1 * mm;
0309   double Envelope_trap_long_edge =
0310       Envelope_trap_short_edge + Envelope_trap_length * tan(prism_angle);
0311
0312   Trap Envelope_trap = MakeTrap("Envelope_trap", Envelope_trap_width, Envelope_trap_length,
0313                                 Envelope_trap_long_edge, Envelope_trap_short_edge);
0314   Position Envelope_trap_position(prism_position_x, 0, prism_position_z - 0.5 * mm);
0315   Volume Envelope_trap_vol("Envelope_trap_vol", Envelope_trap, desc.material("AirOptical"));
0316   dirc_module.placeVolume(Envelope_trap_vol, Transform3D(prism_rotation, Envelope_trap_position));
0317   Envelope_trap_vol.placeVolume(prism_vol, Position(0, 0.5 * mm, 0));
0318
0319   //---- MCP variables
0320   xml_comp_t xml_mcp   = xml_module.child(_Unicode(mcp));
0321   double mcp_thickness = xml_mcp.thickness();
0322   double mcp_height    = xml_mcp.height();
0323   double mcp_width     = xml_mcp.width();
0324   //
0325   //---- MCP construction
0326   Box mcp_box("mcp_box", mcp_height / 2, mcp_width / 2, mcp_thickness / 2);
0327   Volume mcp_vol("mcp_vol", mcp_box, desc.material(xml_mcp.materialStr()));
0328   mcp_vol.setVisAttributes(desc.visAttributes(xml_mcp.visStr())).setSensitiveDetector(sens);
0329   double mcp_position_z = -0.5 * prism_length;
0330   printout(DEBUG, "DIRC_geo", "mcp_position_z LAB   = %12.4f ",
0331            mcp_position_z + dirc_pos_z + Envelope_trap_position.z());
0332   //
0333   Position mcp_position(prism_position_x, mcp_position_z, 0);
0334   RotationX mcp_rotation(-M_PI / 2.);
0335   //
0336   Envelope_trap_vol.placeVolume(mcp_vol, Transform3D(mcp_rotation, mcp_position));
0337
0338   //---- Place modules -----------------------------------------------
0339   const int module_repeat = xml_module.repeat();
0340   const double dphi       = 2. * M_PI / module_repeat;
0341   for (int i = 0; i < module_repeat; i++) {
0342     double phi = dphi * i;
0343     double x   = det_ravg * cos(phi);
0344     double y   = det_ravg * sin(phi);
0345     Transform3D tr(RotationZ(phi), Position(x, y, 0));
0346     det_volume.placeVolume(dirc_module, tr).addPhysVolID("module", i);
0347   }
0348
0349   //---- Construct Assembly dirc_support
0350   //
0351   //    dirc_support contains
0352   //            "rail_trap" -> Volume "rail_vol"
0353   //
0354   xml_comp_t xml_support = xml_det.child(_U(support));
0355   Assembly dirc_support("DIRCSupport");
0356   dirc_support.setVisAttributes(desc.visAttributes(xml_support.visStr()));
0357
0358   //---- Rail
0359   xml_comp_t xml_rail = xml_support.child(_Unicode(rail));
0360   xml_dim_t rail_pos  = xml_rail.position();
0361
0362   double rail_height = (xml_rail.rmax() - xml_rail.rmin()) * cos(dphi / 2);
0363   double rail_width_at_rmax =
0364       xml_rail.width() + 2 * (xml_rail.rmax() - xml_rail.rmin()) * sin(dphi / 2);
0365   double rail_width_at_rmin = xml_rail.width();
0366   double rail_length        = xml_rail.length();
0367   Trap rail_trap("rail_trap", rail_length / 2, 0, 0, rail_height / 2, rail_width_at_rmin / 2,
0368                  rail_width_at_rmax / 2, 0, rail_height / 2, rail_width_at_rmin / 2,
0369                  rail_width_at_rmax / 2, 0);
0370   Volume rail_vol("rail_vol", rail_trap, desc.material(xml_rail.materialStr()));
0371   rail_vol.setVisAttributes(desc.visAttributes(xml_rail.visStr()));
0372
0373   //---- Place rail inside dirc_support
0374   Position rail_position(rail_pos.x(), rail_pos.y(), rail_pos.z());
0375   RotationZ rail_rotation(-M_PI / 2.);
0376   dirc_support.placeVolume(rail_vol, Transform3D(rail_rotation, rail_position));
0377
0378   //---- Place copies of dirc_support in det_volume
0379   for (int i = 0; i < module_repeat; i++) {
0380     double phi = dphi * i + dphi / 2;
0381     double x   = RailRadius * cos(phi);
0382     double y   = RailRadius * sin(phi);
0383     Transform3D tr(RotationZ(phi), Position(x, y, 0));
0384     det_volume.placeVolume(dirc_support, tr);
0385     //
0386     printout(DEBUG, "DIRC_geo", "support. i= %2d \t phi= %6.2f \t x=%12.4f \t y=%12.4f", i, phi, x,
0387              y);
0388     //
0389   }
0390
0391   printout(DEBUG, "DIRC_geo", "---- end hpDIRC geometry ----");
0392
0393 //---- return printlevel to initial value
0394 #ifdef VERBOSE_DIRCgeo
0395   setPrintLevel(priorPrintLevel);
0396 #endif
0397
0398   return det;
0399 }
0400
0401 static dd4hep::Trap MakeTrap(const std::string& pName, double pZ, double pY, double pX,
0402                              double pLTX) {
0403   // Fixed Trap constructor. This function is a workaround of this bug:
0404   // https://github.com/AIDASoft/DD4hep/issues/850
0405   // Should be used instead of dd4hep::Trap(pName, pZ, pY, pX, pLTX) constructor
0406
0407   double fDz         = 0.5 * pZ;
0408   double fTthetaCphi = 0;
0409   double fTthetaSphi = 0;
0410   double fDy1        = 0.5 * pY;
0411   double fDx1        = 0.5 * pX;
0412   double fDx2        = 0.5 * pLTX;
0413   double fTalpha1    = atan(0.5 * (pLTX - pX) / pY);
0414   double fDy2        = fDy1;
0415   double fDx3        = fDx1;
0416   double fDx4        = fDx2;
0417   double fTalpha2    = fTalpha1;
0418
0419   return Trap(pName, fDz, fTthetaCphi, fTthetaSphi, fDy1, fDx1, fDx2, fTalpha1, fDy2, fDx3, fDx4,
0420               fTalpha2);
0421 }
0422
0423 DECLARE_DETELEMENT(epic_DIRC, createDetector)