Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2026-05-07 08:04:39 

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 
0125   // module placement
0126   xml::Component plm = detElem.child(_Unicode(placements));
0127 
0128   std::map<int, std::pair<int, int>> sectorModuleRowsColumns;
0129   auto addRowColumnNumbers = [&sectorModuleRowsColumns](int sector, std::pair<int, int> rowcolumn) {
0130     auto it = sectorModuleRowsColumns.find(sector);
0131     if (it != sectorModuleRowsColumns.end()) {
0132       it->second = rowcolumn;
0133     } else {
0134       sectorModuleRowsColumns[sector] = rowcolumn;
0135     }
0136   };
0137 
0138   int sector_id = 1;
0139   for (xml::Collection_t disk_12surface(plm, _Unicode(disk_12surface)); disk_12surface;
0140        ++disk_12surface) {
0141     auto [sector, rowcolumn] =
0142         add_12surface_disk(desc, assembly, disk_12surface, sens, sector_id++);
0143     addRowColumnNumbers(sector, rowcolumn);
0144   }
0145 
0146   for (auto [sector, rowcolumn] : sectorModuleRowsColumns) {
0147     desc.add(Constant(Form((detName + "_NModules_Sector%d").c_str(), sector),
0148                       std::to_string((rowcolumn.first)), std::to_string((rowcolumn.second))));
0149   }
0150 
0151   // detector position and rotation
0152   auto pos         = get_xml_xyz(detElem, _Unicode(position));
0153   auto rot         = get_xml_xyz(detElem, _Unicode(rotation));
0154   Volume motherVol = desc.pickMotherVolume(det);
0155   Transform3D tr =
0156       Translation3D(pos.x(), pos.y(), pos.z()) * RotationZYX(rot.z(), rot.y(), rot.x());
0157   PlacedVolume envPV = motherVol.placeVolume(assembly, tr);
0158   envPV.addPhysVolID("system", detID);
0159   det.setPlacement(envPV);
0160   return det;
0161 }
0162 
0163 // helper function to build module with or w/o wrapper
0164 std::tuple<Volume, Position> build_module(Detector& desc, xml::Collection_t& plm,
0165                                           SensitiveDetector& sens) {
0166   auto mod = plm.child(_Unicode(module));
0167   auto mx  = mod.attr<double>(_Unicode(modulex));
0168   auto my  = mod.attr<double>(_Unicode(moduley));
0169   auto mz  = mod.attr<double>(_Unicode(modulez));
0170   auto mdz = mod.attr<double>(_Unicode(moduleshift));
0171   Box modshape(mx / 2., my / 2., mz / 2.);
0172   auto modMat = desc.material(mod.attr<std::string>(_Unicode(gmaterial)));
0173   Volume modVol("module_vol", modshape, modMat);
0174   modVol.setVisAttributes(desc.visAttributes(mod.attr<std::string>(_Unicode(vis))));
0175 
0176   auto cry         = plm.child(_Unicode(crystal));
0177   auto cryx        = cry.attr<double>(_Unicode(sizex));
0178   auto cryy        = cry.attr<double>(_Unicode(sizey));
0179   auto cryz        = cry.attr<double>(_Unicode(sizez));
0180   auto roc         = plm.child(_Unicode(readout));
0181   auto PCBx        = roc.attr<double>(_Unicode(PCB_sizex));
0182   auto PCBy        = roc.attr<double>(_Unicode(PCB_sizex));
0183   auto PCBz        = roc.attr<double>(_Unicode(PCB_thickness));
0184   auto sensorx     = roc.attr<double>(_Unicode(Sensor_sizex));
0185   auto sensory     = roc.attr<double>(_Unicode(Sensor_sizey));
0186   auto sensorz     = roc.attr<double>(_Unicode(Sensor_thickness));
0187   auto sensorspace = roc.attr<double>(_Unicode(Sensor_space));
0188   auto sensorNx    = roc.attr<int>(_Unicode(Nsensor_X));
0189   auto sensorNy    = roc.attr<int>(_Unicode(Nsensor_Y));
0190 
0191   Box crystalshape(cryx / 2., cryy / 2., cryz / 2.);
0192   auto crystalMat = desc.material(cry.attr<std::string>(_Unicode(material)));
0193   Volume crystalVol("crystal_vol", crystalshape, crystalMat);
0194   modVol.placeVolume(crystalVol, Position(0., 0., PCBz + sensorz + (cryz - mz) / 2.));
0195   crystalVol.setVisAttributes(desc.visAttributes(cry.attr<std::string>(_Unicode(cryvis))));
0196   crystalVol.setSensitiveDetector(sens);
0197 
0198   Box PCBshape(PCBx / 2., PCBy / 2., PCBz / 2.);
0199   auto PCBMat = desc.material(roc.attr<std::string>(_Unicode(material)));
0200   Volume PCBVol("PCB_vol", PCBshape, PCBMat);
0201   modVol.placeVolume(PCBVol, Position(0., 0., (PCBz - mz) / 2.));
0202 
0203   Box sensorshape(sensorx / 2., sensory / 2., sensorz / 2.);
0204   auto sensorMat = desc.material(roc.attr<std::string>(_Unicode(material)));
0205   Volume sensorVol("sensor_vol", sensorshape, sensorMat);
0206   auto marginx = (PCBx - sensorNx * sensorx - (sensorNx - 1) * sensorspace) / 2.;
0207   auto marginy = (PCBy - sensorNy * sensory - (sensorNy - 1) * sensorspace) / 2.;
0208   auto x0      = marginx + sensorx / 2. - PCBx / 2.;
0209   auto y0      = marginy + sensory / 2. - PCBy / 2.;
0210   for (int i = 0; i < sensorNx; i++)
0211     for (int j = 0; j < sensorNy; j++)
0212       modVol.placeVolume(sensorVol,
0213                          Position(x0 + (sensorx + sensorspace) * i,
0214                                   y0 + (sensory + sensorspace) * j, PCBz + (sensorz - mz) / 2.));
0215 
0216   if (!plm.hasChild(_Unicode(wrapper))) { // no wrapper
0217     printout(DEBUG, "HomogeneousCalorimeter", "without wrapper");
0218     return std::make_tuple(modVol, Position{mx, my, mz});
0219   } else {                                   // build wrapper
0220     auto wrp = plm.child(_Unicode(wrapper)); // Read all the contents in the wrapper block
0221     auto wrapcfthickness = wrp.attr<double>(_Unicode(carbonfiber_thickness));
0222     auto wrapcflength    = wrp.attr<double>(_Unicode(carbonfiber_length));
0223     auto wrapVMthickness = wrp.attr<double>(_Unicode(VM2000_thickness));
0224     auto carbonMat       = desc.material(wrp.attr<std::string>(_Unicode(material_carbon)));
0225     auto wrpMat          = desc.material(wrp.attr<std::string>(_Unicode(material_wrap)));
0226     auto gapMat          = desc.material(wrp.attr<std::string>(_Unicode(material_gap)));
0227 
0228     if (wrapcfthickness < 1e-12 * mm)
0229       return std::make_tuple(modVol, Position{mx, my, mz});
0230 
0231     Box carbonShape(mx / 2., my / 2., wrapcflength / 2.);
0232     Box carbonShape_sub((mx - 2. * wrapcfthickness) / 2., (my - 2. * wrapcfthickness) / 2.,
0233                         wrapcflength / 2.);
0234     SubtractionSolid carbon_subtract(carbonShape, carbonShape_sub, Position(0., 0., 0.));
0235 
0236     Box gapShape(mx / 2., my / 2., (cryz - 2. * wrapcflength) / 2.);
0237     Box gapShape_sub((mx - 2. * wrapcfthickness) / 2., (my - 2. * wrapcfthickness) / 2.,
0238                      (cryz - 2. * wrapcflength) / 2.);
0239     SubtractionSolid gap_subtract(gapShape, gapShape_sub, Position(0., 0., 0.));
0240 
0241     Box wrpVM2000((mx - 2. * wrapcfthickness) / 2., (my - 2. * wrapcfthickness) / 2.,
0242                   (cryz + mdz) / 2.);
0243     Box wrpVM2000_sub((mx - 2. * wrapcfthickness - 2. * wrapVMthickness) / 2.,
0244                       (my - 2. * wrapcfthickness - 2. * wrapVMthickness) / 2., cryz / 2.);
0245     SubtractionSolid wrpVM2000_subtract(wrpVM2000, wrpVM2000_sub, Position(0., 0., -mdz / 2.));
0246 
0247     Volume carbonVol("carbon_vol", carbon_subtract, carbonMat);
0248     Volume gapVol("gap_vol", gap_subtract, gapMat);
0249     Volume wrpVol("wrapper_vol", wrpVM2000_subtract, wrpMat);
0250 
0251     modVol.placeVolume(carbonVol, Position(0., 0.,
0252                                            PCBz + sensorz +
0253                                                (wrapcflength - mz) /
0254                                                    2.)); // put the wrap in the both ends of crystal
0255     modVol.placeVolume(carbonVol,
0256                        Position(0., 0., PCBz + sensorz + cryz - (wrapcflength + mz) / 2.));
0257     modVol.placeVolume(
0258         gapVol,
0259         Position(0., 0.,
0260                  PCBz + sensorz + (cryz - mz) / 2.)); // put the gap between two carbon fiber
0261     modVol.placeVolume(wrpVol, Position(0., 0., PCBz + sensorz + (cryz + mdz - mz) / 2.));
0262 
0263     carbonVol.setVisAttributes(desc.visAttributes(wrp.attr<std::string>(_Unicode(vis_carbon))));
0264     gapVol.setVisAttributes(desc.visAttributes(wrp.attr<std::string>(_Unicode(vis_gap))));
0265     wrpVol.setVisAttributes(desc.visAttributes(wrp.attr<std::string>(_Unicode(vis_wrap))));
0266 
0267     printout(DEBUG, "HomogeneousCalorimeter", "with wrapper");
0268 
0269     return std::make_tuple(modVol, Position{mx, my, mz});
0270   }
0271 }
0272 
0273 // place 12 surface disk of modules
0274 static std::tuple<int, std::pair<int, int>> add_12surface_disk(Detector& desc, Assembly& env,
0275                                                                xml::Collection_t& plm,
0276                                                                SensitiveDetector& sens, int sid) {
0277   auto [modVol, modSize]        = build_module(desc, plm, sens);
0278   int sector_id                 = dd4hep::getAttrOrDefault<int>(plm, _Unicode(sector), sid);
0279   double rmax                   = plm.attr<double>(_Unicode(rmax));
0280   double r12min                 = plm.attr<double>(_Unicode(r12min));
0281   double r12max                 = plm.attr<double>(_Unicode(r12max));
0282   double structure_frame_length = plm.attr<double>(_Unicode(outerringlength));
0283   double calo_module_length     = plm.attr<double>(_Unicode(modulelength));
0284   double envelope_length        = plm.attr<double>(_Unicode(envelope_length));
0285   double Prot                   = plm.attr<double>(_Unicode(protate));
0286   double Nrot                   = plm.attr<double>(_Unicode(nrotate));
0287   double Oring_shift            = plm.attr<double>(_Unicode(outerringshift));
0288   double Innera                 = plm.attr<double>(_Unicode(inneradiusa));
0289   double Innerb                 = plm.attr<double>(_Unicode(inneradiusb));
0290   double phimin                 = dd4hep::getAttrOrDefault<double>(plm, _Unicode(phimin), 0.);
0291   double phimax = dd4hep::getAttrOrDefault<double>(plm, _Unicode(phimax), 2. * M_PI);
0292   xml_coll_t pts_extrudedpolygon(plm, _Unicode(points_extrudedpolygon));
0293 
0294   double half_modx = modSize.x() * 0.5, half_mody = modSize.y() * 0.5;
0295 
0296   //=========================================================
0297   // optional envelope volume and the supporting frame
0298   //=========================================================
0299 
0300   // Material for the structure and mother space
0301   //
0302   Material outer_ring_material =
0303       desc.material(getAttrOrDefault<std::string>(plm, _U(material), "StainlessSteelSAE304"));
0304 
0305   //==============================
0306   // Outer supporting frame
0307   //==============================
0308 
0309   PolyhedraRegular solid_ring12(12, r12min, r12max, structure_frame_length);
0310   Volume ring12_vol("ring12", solid_ring12, outer_ring_material);
0311   Transform3D tr_global_Oring = RotationZYX(Prot, 0., 0.) * Translation3D(0., 0., Oring_shift);
0312   ring12_vol.setVisAttributes(desc.visAttributes(plm.attr<std::string>(_Unicode(vis_struc))));
0313 
0314   //=============================
0315   // The mother volume of modules
0316   //=============================
0317   bool has_envelope = dd4hep::getAttrOrDefault<bool>(plm, _Unicode(envelope), false);
0318   PolyhedraRegular solid_world(12, 0., r12min, envelope_length);
0319   EllipticalTube solid_sub(Innera, Innerb, envelope_length / 2.);
0320   Transform3D subtract_pos = RotationZYX(Nrot, 0., 0.) * Translation3D(1 * cm, 0., 0.);
0321   SubtractionSolid calo_subtract(solid_world, solid_sub, subtract_pos);
0322   Volume env_vol(std::string(env.name()) + "_envelope", calo_subtract, desc.material("Air"));
0323   Transform3D tr_global = RotationZYX(Prot, 0., 0.) * Translation3D(0., 0., 0.);
0324   env_vol.setVisAttributes(desc.visAttributes(plm.attr<std::string>(_Unicode(vis_steel_gap))));
0325 
0326   // Place frames and mother volume of modules into the world volume
0327   //
0328   if (has_envelope) {
0329     env.placeVolume(env_vol, tr_global);          // Place the mother volume for all modules
0330     env.placeVolume(ring12_vol, tr_global_Oring); // Place the outer supporting frame
0331 
0332     xml_comp_t collar_comp   = plm.child(_Unicode(inner_support));
0333     Volume inner_support_vol = build_inner_support(desc, collar_comp, pts_extrudedpolygon);
0334     env_vol.placeVolume(inner_support_vol,
0335                         Transform3D{RotationZ{Nrot}} * Translation3D(collar_comp.x_offset(0.),
0336                                                                      collar_comp.y_offset(0.),
0337                                                                      collar_comp.z_offset(0.)));
0338   }
0339 
0340   //=====================================================================
0341   // Placing The Modules
0342   //=====================================================================
0343 
0344   xml_comp_t placement = plm.child(_Unicode(placement));
0345   auto points =
0346       epic::geo::fillRectangles({placement.x_offset(0.), placement.y_offset(0.)}, modSize.x(),
0347                                 modSize.y(), 0., (rmax / std::cos(Prot)), phimin, phimax);
0348 
0349   std::pair<double, double> c1(0., 0.);
0350   auto polyVertex = epic::geo::getPolygonVertices(c1, (rmax / std::cos(Prot)), M_PI / 12., 12);
0351   std::vector<epic::geo::Point> out_vertices, in_vertices;
0352   for (auto p : polyVertex) {
0353     epic::geo::Point a = {p.first, p.second};
0354     out_vertices.push_back(a);
0355   }
0356 
0357   for (xml_coll_t position_i(pts_extrudedpolygon, _U(position)); position_i; ++position_i) {
0358     xml_comp_t position_comp = position_i;
0359     epic::geo::Point inpt    = {position_comp.x(), position_comp.y()};
0360     in_vertices.push_back(inpt);
0361   }
0362 
0363   double minX = 0., maxX = 0., minY = 0., maxY = 0.;
0364   for (auto& square : points) {
0365     epic::geo::Point box[4] = {{square.x() + half_modx, square.y() + half_mody},
0366                                {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     if (epic::geo::isBoxTotalInsidePolygon(box, out_vertices)) {
0370       if (square.x() < minX)
0371         minX = square.x();
0372       if (square.y() < minY)
0373         minY = square.y();
0374       if (square.x() > maxX)
0375         maxX = square.x();
0376       if (square.y() > maxY)
0377         maxY = square.y();
0378     }
0379   }
0380 
0381   int total_count = 0;
0382   int row = 0, column = 0;
0383   int N_row      = std::round((maxY - minY) / modSize.y());
0384   int N_column   = std::round((maxX - minX) / modSize.x());
0385   auto rowcolumn = std::make_pair(N_row, N_column);
0386 
0387   for (auto& square : points) {
0388     epic::geo::Point box[4] = {{square.x() + half_modx, square.y() + half_mody},
0389                                {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 
0393     if (epic::geo::isBoxTotalInsidePolygon(box, out_vertices)) {
0394       if (!epic::geo::isBoxTotalInsidePolygon(box, in_vertices)) {
0395         column = std::round((square.x() - minX) / modSize.x());
0396         row    = std::round((maxY - square.y()) / modSize.y());
0397         Transform3D tr_local =
0398             RotationZYX(Nrot, 0.0, 0.0) *
0399             Translation3D(square.x(), square.y(),
0400                           std::max((envelope_length - calo_module_length) / 2, 0.));
0401         auto modPV = (has_envelope ? env_vol.placeVolume(modVol, tr_local)
0402                                    : env.placeVolume(modVol, tr_global * tr_local));
0403         modPV.addPhysVolID("sector", sector_id)
0404             .addPhysVolID("row", row)
0405             .addPhysVolID("column", column);
0406         total_count++;
0407       }
0408     }
0409   }
0410 
0411   printout(DEBUG, "HomogeneousCalorimeter_geo", "Number of modules: %d", total_count);
0412   printout(DEBUG, "HomogeneousCalorimeter_geo", "Min X, Y position of module: %.2f, %.2f", minX,
0413            minY);
0414   printout(DEBUG, "HomogeneousCalorimeter_geo", "Max X, Y position of module: %.2f, %.2f", maxX,
0415            maxY);
0416 
0417   return {sector_id, rowcolumn};
0418 }
0419 
0420 //@}
0421 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!