epic/src/DIRC_geo.cpp

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 Wouter Deconinck, Dmitry Romanov
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 static Ref_t createDetector(Detector& desc, xml_h e, SensitiveDetector sens) {
0022   xml_det_t xml_det = e;
0023
0024   // Detector element
0025   string det_name = xml_det.nameStr();
0026   int det_id      = xml_det.id();
0027   DetElement det(det_name, det_id);
0028
0029   // Detector dimension, position, rotation
0030   xml_dim_t dirc_dim = xml_det.dimensions();
0031   xml_dim_t dirc_pos = xml_det.position();
0032   double det_rmin    = dirc_dim.rmin();
0033   double det_rmax    = dirc_dim.rmax();
0034   double det_ravg    = (det_rmin + det_rmax) / 2;
0035
0036   // Detector type
0037   sens.setType("tracker");
0038
0039   // Entire DIRC assembly
0040   Assembly det_volume("DIRC");
0041   det_volume.setVisAttributes(desc.visAttributes(xml_det.visStr()));
0042   Transform3D det_tr(RotationY(0), Position(0.0, 0.0, dirc_pos.z()));
0043   det.setPlacement(
0044       desc.pickMotherVolume(det).placeVolume(det_volume, det_tr).addPhysVolID("system", det_id));
0045
0046   // Construct module
0047   xml_comp_t xml_module = xml_det.child(_U(module));
0048
0049   Assembly dirc_module("DIRCModule");
0050   dirc_module.setVisAttributes(desc.visAttributes(xml_module.visStr()));
0051
0052   // Bar
0053   xml_comp_t xml_bar = xml_module.child(_Unicode(bar));
0054   double bar_height  = xml_bar.height();
0055   double bar_width   = xml_bar.width();
0056   double bar_length  = xml_bar.length();
0057   Box bar_box("bar_box", bar_height / 2, bar_width / 2, bar_length / 2);
0058   Volume bar_vol("bar_vol", bar_box, desc.material(xml_bar.materialStr()));
0059   bar_vol.setVisAttributes(desc.visAttributes(xml_bar.visStr()));
0060
0061   // Glue
0062   xml_comp_t xml_glue   = xml_module.child(_Unicode(glue));
0063   double glue_thickness = xml_glue.thickness();
0064   Box glue_box("glue_box", bar_height / 2, bar_width / 2, glue_thickness / 2);
0065   Volume glue_vol("glue_vol", glue_box, desc.material(xml_glue.materialStr()));
0066   glue_vol.setVisAttributes(desc.visAttributes(xml_glue.visStr()));
0067
0068   auto bar_repeat_y    = xml_bar.attr<int>(_Unicode(repeat_y));
0069   auto bar_repeat_z    = xml_bar.attr<int>(_Unicode(repeat_z));
0070   auto bar_gap         = xml_bar.gap();
0071   auto bar_assm_width  = (bar_width + bar_gap) * bar_repeat_y - bar_gap;
0072   auto bar_assm_length = (bar_length + glue_thickness) * bar_repeat_z;
0073
0074   // Mirror construction
0075   xml_comp_t xml_mirror = xml_module.child(_Unicode(mirror));
0076   auto mirror_width     = xml_mirror.width();
0077   auto mirror_height    = xml_mirror.height();
0078   auto mirror_thickness = xml_mirror.thickness();
0079   Box mirror_box("mirror_box", mirror_height / 2, mirror_width / 2, mirror_thickness / 2);
0080   Volume mirror_vol("mirror_vol", mirror_box, desc.material(xml_mirror.materialStr()));
0081   mirror_vol.setVisAttributes(desc.visAttributes(xml_mirror.visStr()));
0082
0083   // Mirror optical surface
0084   auto surfMgr = desc.surfaceManager();
0085   auto surf    = surfMgr.opticalSurface("DIRC_MirrorOpticalSurface");
0086   SkinSurface skin(desc, det, Form("dirc_mirror_optical_surface"), surf, mirror_vol);
0087   skin.isValid();
0088
0089   // Envelope for bars + mirror
0090   Box Envelope_box("Envelope_box", (mirror_height + 1 * mm) / 2, 5 * (bar_width + 0.15 * mm),
0091                    2 * (bar_length + glue_thickness) + 0.5 * mirror_thickness);
0092   Volume Envelope_box_vol("Envelope_box_vol", Envelope_box, desc.material("AirOptical"));
0093   dirc_module.placeVolume(Envelope_box_vol, Position(0, 0, 0.5 * mirror_thickness));
0094
0095   for (int y_index = 0; y_index < bar_repeat_y; y_index++) {
0096     double y = 0.5 * bar_assm_width - 0.5 * bar_width - (bar_width + bar_gap) * y_index;
0097     for (int z_index = 0; z_index < bar_repeat_z; z_index++) {
0098       double z = 0.5 * bar_assm_length - 0.5 * mirror_thickness - 0.5 * bar_length -
0099                  (bar_length + glue_thickness) * z_index;
0100       Envelope_box_vol.placeVolume(glue_vol,
0101                                    Position(0, y, z - 0.5 * (bar_length + glue_thickness)));
0102       Envelope_box_vol.placeVolume(bar_vol, Position(0, y, z))
0103           .addPhysVolID("section", z_index)
0104           .addPhysVolID("bar", y_index);
0105     }
0106   }
0107
0108   // Place mirror
0109   Envelope_box_vol.placeVolume(mirror_vol, Position(0, 0, 0.5 * bar_assm_length));
0110
0111   // Prism variables
0112   xml_comp_t xml_prism    = xml_module.child(_Unicode(prism));
0113   double prism_angle      = xml_prism.angle();
0114   double prism_width      = xml_prism.width();
0115   double prism_length     = xml_prism.length();
0116   double prism_short_edge = getAttrOrDefault(xml_prism, _Unicode(short_edge), 50 * mm);
0117   double prism_long_edge  = prism_short_edge + prism_length * tan(prism_angle);
0118
0119   // Lens variables
0120   xml_comp_t xml_lens   = xml_module.child(_Unicode(lens));
0121   double lens_shift     = getAttrOrDefault(xml_lens, _Unicode(shift), 0 * mm);
0122   double lens_width     = getAttrOrDefault(xml_lens, _Unicode(width), 35 * mm);
0123   double lens_thickness = getAttrOrDefault(xml_lens, _Unicode(thickness), 12 * mm);
0124   double lens_r1        = getAttrOrDefault(xml_lens, _Unicode(r1), 62 * mm);
0125   double lens_r2        = getAttrOrDefault(xml_lens, _Unicode(r2), 36 * mm);
0126
0127   // Lens construction
0128   // 3-layer spherical lens ---
0129
0130   double lens_radius = sqrt(lens_width * lens_width / 4. + bar_height * bar_height / 4.);
0131
0132   double lens_min_thickness = 2.0 * mm;
0133
0134   double ztrans1 = -lens_thickness / 2. - sqrt(lens_r1 * lens_r1 - lens_radius * lens_radius) +
0135                    lens_min_thickness;
0136   double ztrans2 = -lens_thickness / 2. - sqrt(lens_r2 * lens_r2 - lens_radius * lens_radius) +
0137                    lens_min_thickness * 2;
0138
0139   Box lens_symm_box("lens_symm_box", 0.5 * prism_short_edge, 0.5 * lens_width,
0140                     0.5 * lens_thickness);
0141   Volume Envelope_lens_vol("Envelope_lens_vol", lens_symm_box, desc.material("AirOptical"));
0142
0143   Tube lens_symm_tube(0, lens_radius, 0.5 * lens_thickness);
0144
0145   Sphere lens_sphere1(0, lens_r1);
0146   Sphere lens_sphere2(0, lens_r2);
0147
0148   IntersectionSolid lens_box("lens_box", lens_symm_box, lens_symm_box,
0149                              Position(0, 0, -lens_min_thickness * 2));
0150   IntersectionSolid lens_tube("lens_tube", lens_symm_tube, lens_symm_box,
0151                               Position(0, 0, lens_min_thickness * 2));
0152   UnionSolid lens_box_tube("lens_box_tube", lens_box, lens_tube);
0153
0154   IntersectionSolid lens_layer1_solid("lens_layer1_solid", lens_box_tube, lens_sphere1,
0155                                       Position(0, 0, -ztrans1));
0156   SubtractionSolid lens_layer23_solid("lens_layer23_solid", lens_box_tube, lens_sphere1,
0157                                       Position(0, 0, -ztrans1));
0158
0159   IntersectionSolid lens_layer2_solid("lens_layer2_solid", lens_layer23_solid, lens_sphere2,
0160                                       Position(0, 0, -ztrans2));
0161   SubtractionSolid lens_layer3_solid("lens_layer3_solid", lens_layer23_solid, lens_sphere2,
0162                                      Position(0, 0, -ztrans2));
0163
0164   Volume lens_layer1_vol("lens_layer1_vol", lens_layer1_solid,
0165                          desc.material(xml_lens.attr<std::string>(_Unicode(material1))));
0166   Volume lens_layer2_vol("lens_layer2_vol", lens_layer2_solid,
0167                          desc.material(xml_lens.attr<std::string>(_Unicode(material2))));
0168   Volume lens_layer3_vol("lens_layer3_vol", lens_layer3_solid,
0169                          desc.material(xml_lens.attr<std::string>(_Unicode(material3))));
0170
0171   lens_layer1_vol.setVisAttributes(desc.visAttributes(xml_lens.attr<std::string>(_Unicode(vis1))));
0172   lens_layer2_vol.setVisAttributes(desc.visAttributes(xml_lens.attr<std::string>(_Unicode(vis2))));
0173   lens_layer3_vol.setVisAttributes(desc.visAttributes(xml_lens.attr<std::string>(_Unicode(vis3))));
0174
0175   double lens_position_x = lens_shift;
0176   double lens_position_z = -0.5 * (bar_assm_length + lens_thickness);
0177
0178   for (int y_index = 0; y_index < bar_repeat_y; y_index++) {
0179     double lens_position_y = y_index * lens_width - 0.5 * (prism_width - lens_width);
0180
0181     Position lens_position(lens_position_x, lens_position_y, lens_position_z);
0182     dirc_module.placeVolume(Envelope_lens_vol, lens_position);
0183   }
0184
0185   Envelope_lens_vol.placeVolume(lens_layer1_vol);
0186   Envelope_lens_vol.placeVolume(lens_layer2_vol);
0187   Envelope_lens_vol.placeVolume(lens_layer3_vol);
0188
0189   // Prism construction
0190   Trap prism_trap =
0191       MakeTrap("prism_trap", prism_width, prism_length, prism_long_edge, prism_short_edge);
0192   Volume prism_vol("prism_vol", prism_trap, desc.material(xml_prism.materialStr()));
0193   prism_vol.setVisAttributes(desc.visAttributes(xml_prism.visStr()));
0194
0195   double prism_position_x =
0196       (prism_long_edge + prism_short_edge) / 4. - 0.5 * prism_short_edge + lens_shift;
0197   double prism_position_z = -0.5 * (bar_assm_length + prism_length) - lens_thickness;
0198   RotationX prism_rotation(M_PI / 2.);
0199   Position prism_position(prism_position_x, 0, prism_position_z);
0200
0201   // Envelope for prism + mcp
0202
0203   double Envelope_trap_width      = prism_width + 1 * mm;
0204   double Envelope_trap_length     = prism_length + 1 * mm; // mcp thickness is 1 mm
0205   double Envelope_trap_short_edge = prism_short_edge + 1 * mm;
0206   double Envelope_trap_long_edge =
0207       Envelope_trap_short_edge + Envelope_trap_length * tan(prism_angle);
0208
0209   Trap Envelope_trap = MakeTrap("Envelope_trap", Envelope_trap_width, Envelope_trap_length,
0210                                 Envelope_trap_long_edge, Envelope_trap_short_edge);
0211   Position Envelope_trap_position(prism_position_x, 0, prism_position_z - 0.5 * mm);
0212
0213   Volume Envelope_trap_vol("Envelope_trap_vol", Envelope_trap, desc.material("AirOptical"));
0214   dirc_module.placeVolume(Envelope_trap_vol, Transform3D(prism_rotation, Envelope_trap_position));
0215
0216   Envelope_trap_vol.placeVolume(prism_vol, Position(0, 0.5 * mm, 0));
0217
0218   // MCP variables
0219   xml_comp_t xml_mcp   = xml_module.child(_Unicode(mcp));
0220   double mcp_thickness = xml_mcp.thickness();
0221   double mcp_height    = xml_mcp.height();
0222   double mcp_width     = xml_mcp.width();
0223
0224   // MCP construction
0225   Box mcp_box("mcp_box", mcp_height / 2, mcp_width / 2, mcp_thickness / 2);
0226   Volume mcp_vol("mcp_vol", mcp_box, desc.material(xml_mcp.materialStr()));
0227   mcp_vol.setVisAttributes(desc.visAttributes(xml_mcp.visStr())).setSensitiveDetector(sens);
0228
0229   double mcp_position_z = -0.5 * prism_length;
0230   Position mcp_position(prism_position_x, mcp_position_z, 0);
0231   RotationX mcp_rotation(-M_PI / 2.);
0232
0233   Envelope_trap_vol.placeVolume(mcp_vol, Transform3D(mcp_rotation, mcp_position));
0234
0235   // Place modules
0236   const int module_repeat = xml_module.repeat();
0237   const double dphi       = 2. * M_PI / module_repeat;
0238   for (int i = 0; i < module_repeat; i++) {
0239     double phi = dphi * i;
0240     double x   = det_ravg * cos(phi);
0241     double y   = det_ravg * sin(phi);
0242
0243     Transform3D tr(RotationZ(phi), Position(x, y, 0));
0244     det_volume.placeVolume(dirc_module, tr).addPhysVolID("module", i);
0245   }
0246
0247   // Construct support
0248   xml_comp_t xml_support = xml_det.child(_U(support));
0249   Assembly dirc_support("DIRCSupport");
0250   dirc_support.setVisAttributes(desc.visAttributes(xml_support.visStr()));
0251
0252   // Rail
0253   xml_comp_t xml_rail            = xml_support.child(_Unicode(rail));
0254   xml_dim_t rail_pos             = xml_rail.position();
0255   double rail_height             = xml_rail.height();
0256   double rail_width2             = xml_rail.width();
0257   double rail_distance_to_chord2 = rail_width2 / 2 / tan(dphi / 2);
0258   double rail_distance_to_chord1 = rail_distance_to_chord2 - rail_height;
0259   double rail_width1             = 2 * rail_distance_to_chord1 * tan(dphi / 2);
0260   double rail_length             = xml_rail.length();
0261   Trap rail_trap("rail_trap", rail_length / 2, 0, 0, rail_height / 2, rail_width1 / 2,
0262                  rail_width2 / 2, 0, rail_height / 2, rail_width1 / 2, rail_width2 / 2, 0);
0263   Volume rail_vol("rail_vol", rail_trap, desc.material(xml_rail.materialStr()));
0264   rail_vol.setVisAttributes(desc.visAttributes(xml_rail.visStr()));
0265
0266   // Place rail
0267   Position rail_position(rail_pos.x(), rail_pos.y(), rail_pos.z());
0268   RotationZ rail_rotation(-M_PI / 2.);
0269   dirc_support.placeVolume(rail_vol, Transform3D(rail_rotation, rail_position));
0270
0271   // Place support
0272   for (int i = 0; i < module_repeat; i++) {
0273     double phi = dphi * i + dphi / 2;
0274     double x   = det_ravg * cos(phi);
0275     double y   = det_ravg * sin(phi);
0276
0277     Transform3D tr(RotationZ(phi), Position(x, y, 0));
0278     det_volume.placeVolume(dirc_support, tr);
0279   }
0280
0281   return det;
0282 }
0283
0284 static dd4hep::Trap MakeTrap(const std::string& pName, double pZ, double pY, double pX,
0285                              double pLTX) {
0286   // Fixed Trap constructor. This function is a workaround of this bug:
0287   // https://github.com/AIDASoft/DD4hep/issues/850
0288   // Should be used instead of dd4hep::Trap(pName, pZ, pY, pX, pLTX) constructor
0289
0290   double fDz         = 0.5 * pZ;
0291   double fTthetaCphi = 0;
0292   double fTthetaSphi = 0;
0293   double fDy1        = 0.5 * pY;
0294   double fDx1        = 0.5 * pX;
0295   double fDx2        = 0.5 * pLTX;
0296   double fTalpha1    = atan(0.5 * (pLTX - pX) / pY);
0297   double fDy2        = fDy1;
0298   double fDx3        = fDx1;
0299   double fDx4        = fDx2;
0300   double fTalpha2    = fTalpha1;
0301
0302   return Trap(pName, fDz, fTthetaCphi, fTthetaSphi, fDy1, fDx1, fDx2, fTalpha1, fDy2, fDx3, fDx4,
0303               fTalpha2);
0304 }
0305
0306 DECLARE_DETELEMENT(epic_DIRC, createDetector)