Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2026-07-06 07:53:26 

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 - 2024, Chao Peng, Dmitry Romanov, Pu-Kai Wang, Yuwei Zhu, Dmitry Kalinkin
0003 
0004 #include "DD4hep/DetFactoryHelper.h"
0005 #include "DD4hep/Printout.h"
0006 #include "GeometryHelpers.h"
0007 #include <XML/Helper.h>
0008 #include <XML/Utilities.h>
0009 #include <algorithm>
0010 #include <iostream>
0011 #include <math.h>
0012 #include <tuple>
0013 #include <vector>
0014 
0015 using namespace dd4hep;
0016 
0017 /** \addtogroup calorimeters Calorimeters
0018  */
0019 
0020 /** \addtogroup Homogeneous Calorimeter
0021  * \brief Type: **HomogeneousCalorimeter**.
0022  * \ingroup calorimeters
0023  *
0024  * @{
0025  */
0026 
0027 // headers
0028 static std::tuple<int, std::pair<int, int>> add_12surface_disk(Detector& desc, Assembly& env,
0029                                                                xml::Collection_t& plm,
0030                                                                SensitiveDetector& sens, int id);
0031 
0032 // helper function to get x, y, z if defined in a xml component
0033 template <class XmlComp> Position get_xml_xyz(XmlComp& comp, dd4hep::xml::Strng_t name) {
0034   Position pos(0., 0., 0.);
0035   if (comp.hasChild(name)) {
0036     auto child = comp.child(name);
0037     pos.SetX(dd4hep::getAttrOrDefault<double>(child, _Unicode(x), 0.));
0038     pos.SetY(dd4hep::getAttrOrDefault<double>(child, _Unicode(y), 0.));
0039     pos.SetZ(dd4hep::getAttrOrDefault<double>(child, _Unicode(z), 0.));
0040   }
0041   return pos;
0042 }
0043 
0044 static Volume build_inner_support(Detector& desc, xml_comp_t handle,
0045                                   xml_coll_t pts_extrudedpolygon) {
0046   // This consists of two circular tubes joined by straight sections
0047 
0048   Material inner_ring_material = desc.material(handle.materialStr());
0049 
0050   double electron_r      = handle.attr<double>(_Unicode(electron_r));
0051   double proton_r        = handle.attr<double>(_Unicode(proton_r));
0052   double proton_x_offset = handle.attr<double>(_Unicode(proton_x_offset));
0053   double z_length        = handle.z_length();
0054 
0055   std::vector<double> pt_innerframe_x; //The points information for inner supporting frame
0056   std::vector<double> pt_innerframe_y;
0057   for (xml_coll_t position_i(pts_extrudedpolygon, _U(position)); position_i; ++position_i) {
0058     xml_comp_t position_comp = position_i;
0059     pt_innerframe_x.push_back((position_comp.x()));
0060     pt_innerframe_y.push_back((position_comp.y()));
0061   }
0062 
0063   std::vector<double> sec_z  = {-z_length / 2., z_length / 2.};
0064   std::vector<double> sec_x  = {0., 0.};
0065   std::vector<double> sec_y  = {0., 0.};
0066   std::vector<double> zscale = {1., 1.};
0067 
0068   ExtrudedPolygon inner_support_envelope(pt_innerframe_x, pt_innerframe_y, sec_z, sec_x, sec_y,
0069                                          zscale);
0070 
0071   double straight_section_tilt = acos((electron_r - proton_r) / proton_x_offset);
0072   double straight_section_length =
0073       proton_x_offset + (proton_r - electron_r) * cos(straight_section_tilt);
0074   Position straight_section_offset{
0075       straight_section_length / 2 + cos(straight_section_tilt) * electron_r,
0076       0.,
0077       0.,
0078   };
0079 
0080   // use double length for cutout (any value larger than 1 would work)
0081   double cutout_length = 2 * z_length;
0082 
0083   Tube electron_side{0., electron_r, cutout_length / 2};
0084   Tube proton_side{0., proton_r, cutout_length / 2};
0085   Trd1 electron_proton_straight_section{
0086       electron_r * sin(straight_section_tilt),
0087       proton_r * sin(straight_section_tilt),
0088       cutout_length / 2,
0089       straight_section_length / 2,
0090   };
0091   UnionSolid inner_support_hole{
0092       UnionSolid{
0093           electron_side,
0094           proton_side,
0095           Position{proton_x_offset, 0., 0.},
0096       },
0097       electron_proton_straight_section,
0098       Transform3D{straight_section_offset} * RotationZ(90 * deg) * RotationX(90 * deg),
0099   };
0100 
0101   SubtractionSolid inner_support{
0102       inner_support_envelope,
0103       inner_support_hole,
0104   };
0105 
0106   Volume inner_support_vol{"inner_support_vol", inner_support, inner_ring_material};
0107   inner_support_vol.setVisAttributes(desc.visAttributes(handle.visStr()));
0108   return inner_support_vol;
0109 }
0110 
0111 // main
0112 static Ref_t create_detector(Detector& desc, xml::Handle_t handle, SensitiveDetector sens) {
0113   xml::DetElement detElem = handle;
0114   std::string detName     = detElem.nameStr();
0115   int detID               = detElem.id();
0116   DetElement det(detName, detID);
0117   sens.setType("calorimeter");
0118 
0119   // apply any detector type flags set in XML
0120   dd4hep::xml::setDetectorTypeFlag(detElem, det);
0121 
0122   // assembly
0123   Assembly assembly(detName);
0124   assembly.setAttributes(desc, detElem.regionStr(), detElem.limitsStr(), detElem.visStr());
0125 
0126   // module placement
0127   xml::Component plm = detElem.child(_Unicode(placements));
0128 
0129   std::map<int, std::pair<int, int>> sectorModuleRowsColumns;
0130   auto addRowColumnNumbers = [&sectorModuleRowsColumns](int sector, std::pair<int, int> rowcolumn) {
0131     auto it = sectorModuleRowsColumns.find(sector);
0132     if (it != sectorModuleRowsColumns.end()) {
0133       it->second = rowcolumn;
0134     } else {
0135       sectorModuleRowsColumns[sector] = rowcolumn;
0136     }
0137   };
0138 
0139   int sector_id = 1;
0140   for (xml::Collection_t disk_12surface(plm, _Unicode(disk_12surface)); disk_12surface;
0141        ++disk_12surface) {
0142     auto [sector, rowcolumn] =
0143         add_12surface_disk(desc, assembly, disk_12surface, sens, sector_id++);
0144     addRowColumnNumbers(sector, rowcolumn);
0145   }
0146 
0147   for (auto [sector, rowcolumn] : sectorModuleRowsColumns) {
0148     desc.add(Constant(Form((detName + "_NModules_Sector%d").c_str(), sector),
0149                       std::to_string((rowcolumn.first)), std::to_string((rowcolumn.second))));
0150   }
0151 
0152   // detector position and rotation
0153   auto pos         = get_xml_xyz(detElem, _Unicode(position));
0154   auto rot         = get_xml_xyz(detElem, _Unicode(rotation));
0155   Volume motherVol = desc.pickMotherVolume(det);
0156   Transform3D tr =
0157       Translation3D(pos.x(), pos.y(), pos.z()) * RotationZYX(rot.z(), rot.y(), rot.x());
0158   PlacedVolume envPV = motherVol.placeVolume(assembly, tr);
0159   envPV.addPhysVolID("system", detID);
0160   det.setPlacement(envPV);
0161   return det;
0162 }
0163 
0164 // helper function to build module with or w/o wrapper
0165 std::tuple<Volume, Position> build_module(Detector& desc, xml::Collection_t& plm,
0166                                           SensitiveDetector& sens) {
0167   auto mod = plm.child(_Unicode(module));
0168   auto mx  = mod.attr<double>(_Unicode(modulex));
0169   auto my  = mod.attr<double>(_Unicode(moduley));
0170   auto mz  = mod.attr<double>(_Unicode(modulez));
0171   auto mdz = mod.attr<double>(_Unicode(moduleshift));
0172   Box modshape(mx / 2., my / 2., mz / 2.);
0173   auto modMat = desc.material(mod.attr<std::string>(_Unicode(gmaterial)));
0174   Volume modVol("module_vol", modshape, modMat);
0175   modVol.setVisAttributes(desc.visAttributes(mod.attr<std::string>(_Unicode(vis))));
0176 
0177   auto cry         = plm.child(_Unicode(crystal));
0178   auto cryx        = cry.attr<double>(_Unicode(sizex));
0179   auto cryy        = cry.attr<double>(_Unicode(sizey));
0180   auto cryz        = cry.attr<double>(_Unicode(sizez));
0181   auto roc         = plm.child(_Unicode(readout));
0182   auto PCBx        = roc.attr<double>(_Unicode(PCB_sizex));
0183   auto PCBy        = roc.attr<double>(_Unicode(PCB_sizex));
0184   auto PCBz        = roc.attr<double>(_Unicode(PCB_thickness));
0185   auto sensorx     = roc.attr<double>(_Unicode(Sensor_sizex));
0186   auto sensory     = roc.attr<double>(_Unicode(Sensor_sizey));
0187   auto sensorz     = roc.attr<double>(_Unicode(Sensor_thickness));
0188   auto sensorspace = roc.attr<double>(_Unicode(Sensor_space));
0189   auto sensorNx    = roc.attr<int>(_Unicode(Nsensor_X));
0190   auto sensorNy    = roc.attr<int>(_Unicode(Nsensor_Y));
0191 
0192   Box crystalshape(cryx / 2., cryy / 2., cryz / 2.);
0193   auto crystalMat = desc.material(cry.attr<std::string>(_Unicode(material)));
0194   Volume crystalVol("crystal_vol", crystalshape, crystalMat);
0195   modVol.placeVolume(crystalVol, Position(0., 0., PCBz + sensorz + (cryz - mz) / 2.));
0196   crystalVol.setVisAttributes(desc.visAttributes(cry.attr<std::string>(_Unicode(cryvis))));
0197   crystalVol.setSensitiveDetector(sens);
0198 
0199   Box PCBshape(PCBx / 2., PCBy / 2., PCBz / 2.);
0200   auto PCBMat = desc.material(roc.attr<std::string>(_Unicode(material)));
0201   Volume PCBVol("PCB_vol", PCBshape, PCBMat);
0202   modVol.placeVolume(PCBVol, Position(0., 0., (PCBz - mz) / 2.));
0203 
0204   Box sensorshape(sensorx / 2., sensory / 2., sensorz / 2.);
0205   auto sensorMat = desc.material(roc.attr<std::string>(_Unicode(material)));
0206   Volume sensorVol("sensor_vol", sensorshape, sensorMat);
0207   auto marginx = (PCBx - sensorNx * sensorx - (sensorNx - 1) * sensorspace) / 2.;
0208   auto marginy = (PCBy - sensorNy * sensory - (sensorNy - 1) * sensorspace) / 2.;
0209   auto x0      = marginx + sensorx / 2. - PCBx / 2.;
0210   auto y0      = marginy + sensory / 2. - PCBy / 2.;
0211   for (int i = 0; i < sensorNx; i++)
0212     for (int j = 0; j < sensorNy; j++)
0213       modVol.placeVolume(sensorVol,
0214                          Position(x0 + (sensorx + sensorspace) * i,
0215                                   y0 + (sensory + sensorspace) * j, PCBz + (sensorz - mz) / 2.));
0216 
0217   if (!plm.hasChild(_Unicode(wrapper))) { // no wrapper
0218     printout(DEBUG, "HomogeneousCalorimeter", "without wrapper");
0219     return std::make_tuple(modVol, Position{mx, my, mz});
0220   } else {                                   // build wrapper
0221     auto wrp = plm.child(_Unicode(wrapper)); // Read all the contents in the wrapper block
0222     auto wrapcfthickness = wrp.attr<double>(_Unicode(carbonfiber_thickness));
0223     auto wrapcflength    = wrp.attr<double>(_Unicode(carbonfiber_length));
0224     auto wrapVMthickness = wrp.attr<double>(_Unicode(VM2000_thickness));
0225     auto carbonMat       = desc.material(wrp.attr<std::string>(_Unicode(material_carbon)));
0226     auto wrpMat          = desc.material(wrp.attr<std::string>(_Unicode(material_wrap)));
0227     auto gapMat          = desc.material(wrp.attr<std::string>(_Unicode(material_gap)));
0228 
0229     if (wrapcfthickness < 1e-12 * mm)
0230       return std::make_tuple(modVol, Position{mx, my, mz});
0231 
0232     Box carbonShape(mx / 2., my / 2., wrapcflength / 2.);
0233     Box carbonShape_sub((mx - 2. * wrapcfthickness) / 2., (my - 2. * wrapcfthickness) / 2.,
0234                         wrapcflength / 2.);
0235     SubtractionSolid carbon_subtract(carbonShape, carbonShape_sub, Position(0., 0., 0.));
0236 
0237     Box gapShape(mx / 2., my / 2., (cryz - 2. * wrapcflength) / 2.);
0238     Box gapShape_sub((mx - 2. * wrapcfthickness) / 2., (my - 2. * wrapcfthickness) / 2.,
0239                      (cryz - 2. * wrapcflength) / 2.);
0240     SubtractionSolid gap_subtract(gapShape, gapShape_sub, Position(0., 0., 0.));
0241 
0242     Box wrpVM2000((mx - 2. * wrapcfthickness) / 2., (my - 2. * wrapcfthickness) / 2.,
0243                   (cryz + mdz) / 2.);
0244     Box wrpVM2000_sub((mx - 2. * wrapcfthickness - 2. * wrapVMthickness) / 2.,
0245                       (my - 2. * wrapcfthickness - 2. * wrapVMthickness) / 2., cryz / 2.);
0246     SubtractionSolid wrpVM2000_subtract(wrpVM2000, wrpVM2000_sub, Position(0., 0., -mdz / 2.));
0247 
0248     Volume carbonVol("carbon_vol", carbon_subtract, carbonMat);
0249     Volume gapVol("gap_vol", gap_subtract, gapMat);
0250     Volume wrpVol("wrapper_vol", wrpVM2000_subtract, wrpMat);
0251 
0252     modVol.placeVolume(carbonVol, Position(0., 0.,
0253                                            PCBz + sensorz +
0254                                                (wrapcflength - mz) /
0255                                                    2.)); // put the wrap in the both ends of crystal
0256     modVol.placeVolume(carbonVol,
0257                        Position(0., 0., PCBz + sensorz + cryz - (wrapcflength + mz) / 2.));
0258     modVol.placeVolume(
0259         gapVol,
0260         Position(0., 0.,
0261                  PCBz + sensorz + (cryz - mz) / 2.)); // put the gap between two carbon fiber
0262     modVol.placeVolume(wrpVol, Position(0., 0., PCBz + sensorz + (cryz + mdz - mz) / 2.));
0263 
0264     carbonVol.setVisAttributes(desc.visAttributes(wrp.attr<std::string>(_Unicode(vis_carbon))));
0265     gapVol.setVisAttributes(desc.visAttributes(wrp.attr<std::string>(_Unicode(vis_gap))));
0266     wrpVol.setVisAttributes(desc.visAttributes(wrp.attr<std::string>(_Unicode(vis_wrap))));
0267 
0268     printout(DEBUG, "HomogeneousCalorimeter", "with wrapper");
0269 
0270     return std::make_tuple(modVol, Position{mx, my, mz});
0271   }
0272 }
0273 
0274 // place 12 surface disk of modules
0275 static std::tuple<int, std::pair<int, int>> add_12surface_disk(Detector& desc, Assembly& env,
0276                                                                xml::Collection_t& plm,
0277                                                                SensitiveDetector& sens, int sid) {
0278   auto [modVol, modSize]        = build_module(desc, plm, sens);
0279   int sector_id                 = dd4hep::getAttrOrDefault<int>(plm, _Unicode(sector), sid);
0280   double rmax                   = plm.attr<double>(_Unicode(rmax));
0281   double r12min                 = plm.attr<double>(_Unicode(r12min));
0282   double r12max                 = plm.attr<double>(_Unicode(r12max));
0283   double structure_frame_length = plm.attr<double>(_Unicode(outerringlength));
0284   double calo_module_length     = plm.attr<double>(_Unicode(modulelength));
0285   double envelope_length        = plm.attr<double>(_Unicode(envelope_length));
0286   double Prot                   = plm.attr<double>(_Unicode(protate));
0287   double Nrot                   = plm.attr<double>(_Unicode(nrotate));
0288   double Oring_shift            = plm.attr<double>(_Unicode(outerringshift));
0289   double Innera                 = plm.attr<double>(_Unicode(inneradiusa));
0290   double Innerb                 = plm.attr<double>(_Unicode(inneradiusb));
0291   double phimin                 = dd4hep::getAttrOrDefault<double>(plm, _Unicode(phimin), 0.);
0292   double phimax = dd4hep::getAttrOrDefault<double>(plm, _Unicode(phimax), 2. * M_PI);
0293   xml_coll_t pts_extrudedpolygon(plm, _Unicode(points_extrudedpolygon));
0294 
0295   double half_modx = modSize.x() * 0.5, half_mody = modSize.y() * 0.5;
0296 
0297   //=========================================================
0298   // optional envelope volume and the supporting frame
0299   //=========================================================
0300 
0301   // Material for the structure and mother space
0302   //
0303   Material outer_ring_material =
0304       desc.material(getAttrOrDefault<std::string>(plm, _U(material), "StainlessSteelSAE304"));
0305 
0306   //==============================
0307   // Outer supporting frame
0308   //==============================
0309 
0310   PolyhedraRegular solid_ring12(12, r12min, r12max, structure_frame_length);
0311   Volume ring12_vol("ring12", solid_ring12, outer_ring_material);
0312   Transform3D tr_global_Oring = RotationZYX(Prot, 0., 0.) * Translation3D(0., 0., Oring_shift);
0313   ring12_vol.setVisAttributes(desc.visAttributes(plm.attr<std::string>(_Unicode(vis_struc))));
0314 
0315   //=============================
0316   // The mother volume of modules
0317   //=============================
0318   bool has_envelope = dd4hep::getAttrOrDefault<bool>(plm, _Unicode(envelope), false);
0319   PolyhedraRegular solid_world(12, 0., r12min, envelope_length);
0320   EllipticalTube solid_sub(Innera, Innerb, envelope_length / 2.);
0321   Transform3D subtract_pos = RotationZYX(Nrot, 0., 0.) * Translation3D(1 * cm, 0., 0.);
0322   SubtractionSolid calo_subtract(solid_world, solid_sub, subtract_pos);
0323   Volume env_vol(std::string(env.name()) + "_envelope", calo_subtract, desc.material("Air"));
0324   Transform3D tr_global = RotationZYX(Prot, 0., 0.) * Translation3D(0., 0., 0.);
0325   env_vol.setVisAttributes(desc.visAttributes(plm.attr<std::string>(_Unicode(vis_steel_gap))));
0326 
0327   // Place frames and mother volume of modules into the world volume
0328   //
0329   if (has_envelope) {
0330     env.placeVolume(env_vol, tr_global);          // Place the mother volume for all modules
0331     env.placeVolume(ring12_vol, tr_global_Oring); // Place the outer supporting frame
0332 
0333     xml_comp_t collar_comp   = plm.child(_Unicode(inner_support));
0334     Volume inner_support_vol = build_inner_support(desc, collar_comp, pts_extrudedpolygon);
0335     env_vol.placeVolume(inner_support_vol,
0336                         Transform3D{RotationZ{Nrot}} * Translation3D(collar_comp.x_offset(0.),
0337                                                                      collar_comp.y_offset(0.),
0338                                                                      collar_comp.z_offset(0.)));
0339   }
0340 
0341   //=====================================================================
0342   // Placing The Modules
0343   //=====================================================================
0344 
0345   xml_comp_t placement = plm.child(_Unicode(placement));
0346   auto points =
0347       epic::geo::fillRectangles({placement.x_offset(0.), placement.y_offset(0.)}, modSize.x(),
0348                                 modSize.y(), 0., (rmax / std::cos(Prot)), phimin, phimax);
0349 
0350   std::pair<double, double> c1(0., 0.);
0351   auto polyVertex = epic::geo::getPolygonVertices(c1, (rmax / std::cos(Prot)), M_PI / 12., 12);
0352   std::vector<epic::geo::Point> out_vertices, in_vertices;
0353   for (auto p : polyVertex) {
0354     epic::geo::Point a = {p.first, p.second};
0355     out_vertices.push_back(a);
0356   }
0357 
0358   for (xml_coll_t position_i(pts_extrudedpolygon, _U(position)); position_i; ++position_i) {
0359     xml_comp_t position_comp = position_i;
0360     epic::geo::Point inpt    = {position_comp.x(), position_comp.y()};
0361     in_vertices.push_back(inpt);
0362   }
0363 
0364   double minX = 0., maxX = 0., minY = 0., maxY = 0.;
0365   for (auto& square : points) {
0366     epic::geo::Point box[4] = {{square.x() + half_modx, square.y() + half_mody},
0367                                {square.x() - half_modx, square.y() + half_mody},
0368                                {square.x() - half_modx, square.y() - half_mody},
0369                                {square.x() + half_modx, square.y() - half_mody}};
0370     if (epic::geo::isBoxTotalInsidePolygon(box, out_vertices)) {
0371       if (square.x() < minX)
0372         minX = square.x();
0373       if (square.y() < minY)
0374         minY = square.y();
0375       if (square.x() > maxX)
0376         maxX = square.x();
0377       if (square.y() > maxY)
0378         maxY = square.y();
0379     }
0380   }
0381 
0382   int total_count = 0;
0383   int row = 0, column = 0;
0384   int N_row      = std::round((maxY - minY) / modSize.y());
0385   int N_column   = std::round((maxX - minX) / modSize.x());
0386   auto rowcolumn = std::make_pair(N_row, N_column);
0387 
0388   for (auto& square : points) {
0389     epic::geo::Point box[4] = {{square.x() + half_modx, square.y() + half_mody},
0390                                {square.x() - half_modx, square.y() + half_mody},
0391                                {square.x() - half_modx, square.y() - half_mody},
0392                                {square.x() + half_modx, square.y() - half_mody}};
0393 
0394     if (epic::geo::isBoxTotalInsidePolygon(box, out_vertices)) {
0395       if (!epic::geo::isBoxTotalInsidePolygon(box, in_vertices)) {
0396         column = std::round((square.x() - minX) / modSize.x());
0397         row    = std::round((maxY - square.y()) / modSize.y());
0398         Transform3D tr_local =
0399             RotationZYX(Nrot, 0.0, 0.0) *
0400             Translation3D(square.x(), square.y(),
0401                           std::max((envelope_length - calo_module_length) / 2, 0.));
0402         auto modPV = (has_envelope ? env_vol.placeVolume(modVol, tr_local)
0403                                    : env.placeVolume(modVol, tr_global * tr_local));
0404         modPV.addPhysVolID("sector", sector_id)
0405             .addPhysVolID("row", row)
0406             .addPhysVolID("column", column);
0407         total_count++;
0408       }
0409     }
0410   }
0411 
0412   printout(DEBUG, "HomogeneousCalorimeter_geo", "Number of modules: %d", total_count);
0413   printout(DEBUG, "HomogeneousCalorimeter_geo", "Min X, Y position of module: %.2f, %.2f", minX,
0414            minY);
0415   printout(DEBUG, "HomogeneousCalorimeter_geo", "Max X, Y position of module: %.2f, %.2f", maxX,
0416            maxY);
0417 
0418   return {sector_id, rowcolumn};
0419 }
0420 
0421 //@}
0422 DECLARE_DETELEMENT(epic_HomogeneousCalorimeter, create_detector)
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!