epic/src/InsertCalorimeter_geo.cpp

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 Ryan Milton
0003
0004 //==========================================================================
0005 //  Implementation of forward insert calorimeter
0006 //--------------------------------------------------------------------------
0007 //  Author: Ryan Milton (UCR)
0008 //==========================================================================
0009
0010 #include "DD4hep/DetFactoryHelper.h"
0011 #include <XML/Helper.h>
0012 #include <XML/Layering.h>
0013 #include <tuple>
0014 #include <vector>
0015
0016 using namespace dd4hep;
0017
0018 static Ref_t createDetector(Detector& desc, xml_h handle, SensitiveDetector sens) {
0019   xml_det_t detElem   = handle;
0020   std::string detName = detElem.nameStr();
0021   int detID           = detElem.id();
0022
0023   xml_dim_t dim = detElem.dimensions();
0024   double width  = dim.x(); // Size along x-axis
0025   double height = dim.y(); // Size along y-axis
0026   double length = dim.z(); // Size along z-axis
0027
0028   xml_dim_t pos = detElem.position(); // Position in global coordinates
0029
0030   Material air = desc.material("Air");
0031
0032   // Getting beampipe hole dimensions
0033   const xml::Component& beampipe_hole_xml = detElem.child(_Unicode(beampipe_hole));
0034   const double hole_radius_initial        = dd4hep::getAttrOrDefault<double>(
0035       beampipe_hole_xml, _Unicode(initial_hole_radius), 14.61 * cm);
0036   const double hole_radius_final =
0037       dd4hep::getAttrOrDefault<double>(beampipe_hole_xml, _Unicode(final_hole_radius), 17.17 * cm);
0038   const std::pair<double, double> hole_radii_parameters(hole_radius_initial, hole_radius_final);
0039
0040   // Subtract by pos.x() and pos.y() to convert from global to local coordinates
0041   const double hole_x_initial =
0042       dd4hep::getAttrOrDefault<double>(beampipe_hole_xml, _Unicode(initial_hole_x), -7.20 * cm) -
0043       pos.x();
0044   const double hole_x_final =
0045       dd4hep::getAttrOrDefault<double>(beampipe_hole_xml, _Unicode(final_hole_x), -10.44 * cm) -
0046       pos.x();
0047   const std::pair<double, double> hole_x_parameters(hole_x_initial, hole_x_final);
0048
0049   const double hole_y_initial =
0050       dd4hep::getAttrOrDefault<double>(beampipe_hole_xml, _Unicode(initial_hole_y), 0. * cm) -
0051       pos.y();
0052   const double hole_y_final =
0053       dd4hep::getAttrOrDefault<double>(beampipe_hole_xml, _Unicode(final_hole_y), 0. * cm) -
0054       pos.y();
0055   const std::pair<double, double> hole_y_parameters(hole_y_initial, hole_y_final);
0056
0057   const double left_right_gap =
0058       dd4hep::getAttrOrDefault<double>(detElem, _Unicode(left_right_gap), 0.38 * cm);
0059   // Getting thickness of backplate
0060   /*
0061     The hole radius & position is determined by liner interpolation
0062     The HCal insert has multiple layers; interpolation goes from front of insert to front of final layer
0063     So need final layer thickness (backplate thickness) for interpolation
0064
0065     For the ECal insert, the hole radius & position is constant
0066     Also has only one layer so don't have a backplate_thickness there (so set to 0)
0067   */
0068   auto backplate_thickness =
0069       detElem.hasChild(_Unicode(backplate))
0070           ? detElem.child(_Unicode(backplate)).attr<double>(_Unicode(thickness))
0071           : 0.;
0072
0073   // Function that returns a linearly interpolated hole radius, x-position, and y-position at a given z
0074   auto get_hole_rxy = [hole_radii_parameters, hole_x_parameters, hole_y_parameters, length,
0075                        backplate_thickness](double z_pos) {
0076     /*
0077       radius = (hole_radius_final - hole_radius_initial)/(length) * z_pos +  hole_radius_initial
0078       Treats the beginning of the insert as z = 0
0079       At z = 0 (beginning of first layer), hole radius is hole_radius_initial
0080       The radius is hole_radius_final at the beginning of the backplate,
0081         i.e. z = length - backplate_thickness
0082     */
0083     double hole_radius_slope = (hole_radii_parameters.second - hole_radii_parameters.first) /
0084                                (length - backplate_thickness);
0085     double hole_radius_at_z = hole_radius_slope * z_pos + hole_radii_parameters.first;
0086
0087     double hole_xpos_slope =
0088         (hole_x_parameters.second - hole_x_parameters.first) / (length - backplate_thickness);
0089     double hole_xpos_at_z = hole_xpos_slope * z_pos + hole_x_parameters.first;
0090
0091     double hole_ypos_slope =
0092         (hole_y_parameters.second - hole_y_parameters.first) / (length - backplate_thickness);
0093     double hole_ypos_at_z = hole_ypos_slope * z_pos + hole_y_parameters.first;
0094
0095     std::tuple<double, double, double> hole_rxy =
0096         std::make_tuple(hole_radius_at_z, hole_xpos_at_z, hole_ypos_at_z);
0097     return hole_rxy;
0098   };
0099
0100   // Assembly that will contain all the layers
0101   Assembly assembly(detName);
0102   // FIXME Workaround for https://github.com/eic/epic/issues/411
0103   assembly.setVisAttributes(desc.visAttributes("InvisibleWithDaughters"));
0104   PlacedVolume pv;
0105   for (int side_num = 0; side_num < 2; side_num++) { // 0 = right, 1 = left
0106     std::string side_name = side_num == 1 ? "L" : "R";
0107     // Keeps track of the z location as we move longiduinally through the insert
0108     // Will use this tracking variable as input to get_hole_rxy
0109     double z_distance_traversed = 0.;
0110
0111     int layer_num = 1;
0112
0113     // Looping through all the different layer sections (W/Sc, Steel/Sc, backplate)
0114     for (xml_coll_t c(detElem, _U(layer)); c; c++) {
0115       xml_comp_t x_layer     = c;
0116       int repeat             = x_layer.repeat();
0117       double layer_thickness = x_layer.thickness();
0118
0119       // Looping through the number of repeated layers in each section
0120       for (int i = 0; i < repeat; i++) {
0121         std::string layer_name = detName + _toString(layer_num, "_layer%d") + "_" + side_name;
0122         Box layer(width / 2., height / 2., layer_thickness / 2.);
0123
0124         // Hole radius and position for each layer is determined from z position at the front of the layer
0125         const auto hole_rxy = get_hole_rxy(z_distance_traversed);
0126         double hole_r       = std::get<0>(hole_rxy);
0127         double hole_x       = std::get<1>(hole_rxy);
0128         double hole_y       = std::get<2>(hole_rxy);
0129
0130         // Removing beampipe shape from each layer
0131         Tube layer_hole(0., hole_r, layer_thickness / 2.);
0132         SubtractionSolid layer_with_hole(layer, layer_hole, Position(hole_x, hole_y, 0.));
0133         // Only select the left or right side of the layer
0134         Box side_cut(width / 2., height, layer_thickness);
0135         Position side_cut_position((width / 2 - left_right_gap / 2) * (1 - 2 * side_num) - pos.x(),
0136                                    0, 0);
0137         SubtractionSolid layer_side_with_hole(layer_with_hole, side_cut, side_cut_position);
0138         Volume layer_vol(layer_name, layer_side_with_hole, air);
0139
0140         int slice_num  = 1;
0141         double slice_z = -layer_thickness / 2.; // Keeps track of slices' z locations in each layer
0142
0143         // Looping over each layer's slices
0144         for (xml_coll_t l(x_layer, _U(slice)); l; l++) {
0145           xml_comp_t x_slice     = l;
0146           double slice_thickness = x_slice.thickness();
0147           std::string slice_name = layer_name + _toString(slice_num, "slice%d");
0148           Material slice_mat     = desc.material(x_slice.materialStr());
0149           slice_z += slice_thickness / 2.; // Going to slice halfway point
0150
0151           // Each slice within a layer has the same hole radius and x-y position
0152           Box slice(width / 2., height / 2., slice_thickness / 2.);
0153           Tube slice_hole(0., hole_r, slice_thickness / 2.);
0154           SubtractionSolid slice_with_hole(slice, slice_hole, Position(hole_x, hole_y, 0.));
0155           Box side_cut_slice(width / 2., height, layer_thickness);
0156           Position side_cut_position_slice(
0157               (width / 2 - left_right_gap / 2) * (1 - 2 * side_num) - pos.x(), 0, 0);
0158           SubtractionSolid slice_side_with_hole(slice_with_hole, side_cut_slice,
0159                                                 side_cut_position_slice);
0160           Volume slice_vol(slice_name, slice_side_with_hole, slice_mat);
0161
0162           // Setting appropriate slices as sensitive
0163           if (x_slice.isSensitive()) {
0164             sens.setType("calorimeter");
0165             slice_vol.setSensitiveDetector(sens);
0166           }
0167
0168           // Setting slice attributes
0169           slice_vol.setAttributes(desc, x_slice.regionStr(), x_slice.limitsStr(), x_slice.visStr());
0170
0171           // Placing slice within layer
0172           pv = layer_vol.placeVolume(slice_vol,
0173                                      Transform3D(RotationZYX(0, 0, 0), Position(0., 0., slice_z)));
0174           pv.addPhysVolID("slice", slice_num);
0175           pv.addPhysVolID("side", side_num);
0176           slice_z += slice_thickness / 2.;
0177           z_distance_traversed += slice_thickness;
0178           slice_num++;
0179         }
0180
0181         // Setting layer attributes
0182         layer_vol.setAttributes(desc, x_layer.regionStr(), x_layer.limitsStr(), x_layer.visStr());
0183         /*
0184           Placing each layer inside assembly
0185           -length/2. is front of detector in global coordinate system
0186           + (z_distance_traversed - layer_thickness) goes to the front of each layer
0187           + layer_thickness/2. places layer in correct spot
0188           Example: After placement of slices in first layer, z_distance_traversed = layer_thickness
0189                    Subtracting layer_thickness goes back to the front of the first slice (Now, z = -length/2)
0190                    Adding layer_thickness / 2. goes to half the first layer thickness (proper place to put layer)
0191                    Each loop over repeat will increases z_distance_traversed by layer_thickness
0192         */
0193         pv = assembly.placeVolume(
0194             layer_vol,
0195             Transform3D(RotationZYX(0, 0, 0),
0196                         Position(0., 0.,
0197                                  -length / 2. + (z_distance_traversed - layer_thickness) +
0198                                      layer_thickness / 2.)));
0199
0200         pv.addPhysVolID("layer", layer_num);
0201         pv.addPhysVolID("side", side_num);
0202         layer_num++;
0203       }
0204     }
0205   }
0206   DetElement det(detName, detID);
0207   Volume motherVol = desc.pickMotherVolume(det);
0208
0209   // Placing insert in world volume
0210   auto tr          = Transform3D(Position(pos.x(), pos.y(), pos.z() + length / 2.));
0211   PlacedVolume phv = motherVol.placeVolume(assembly, tr);
0212   phv.addPhysVolID("system", detID);
0213   det.setPlacement(phv);
0214
0215   return det;
0216 }
0217 DECLARE_DETELEMENT(epic_InsertCalorimeter, createDetector)