Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-01-18 09:15:57 

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