epic/src/CryostatMagnet_geo.cpp

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 Alex Jentsch, Wouter Deconinck, Whitney Armstrong, Andrii Natochii
0003
0004 #include "DD4hep/DetFactoryHelper.h"
0005 #include "DD4hep/Printout.h"
0006 #include "DD4hep/Shapes.h"
0007 #include "DDRec/DetectorData.h"
0008 #include "DDRec/Surface.h"
0009 #include "TMath.h"
0010 #include "XML/Layering.h"
0011
0012 using namespace std;
0013 using namespace dd4hep;
0014 using namespace dd4hep::rec;
0015 using namespace ROOT::Math;
0016
0017 void buildTubeElement(dd4hep::DetElement& sdet, dd4hep::Assembly& assembly, dd4hep::Detector& dtor,
0018                       dd4hep::xml::DetElement x_det, string element, dd4hep::Material material);
0019 void buildPolyElement(dd4hep::DetElement& sdet, dd4hep::Assembly& assembly, dd4hep::Detector& dtor,
0020                       dd4hep::xml::DetElement x_det, string element, dd4hep::Material material);
0021
0022 static Ref_t build_magnet(Detector& dtor, xml_h e, SensitiveDetector /* sens */) {
0023   xml_det_t x_det = e;
0024   int det_id      = x_det.id();
0025   string det_name = x_det.nameStr();
0026
0027   DetElement sdet(det_name, det_id);
0028   Assembly assembly(det_name + "_assembly");
0029
0030   // get materials
0031   Material iron       = dtor.material("Iron");
0032   Material nbti       = dtor.material("SolenoidCoil");
0033   Material steel_304l = dtor.material("StainlessSteelSAE304");
0034   Material alum       = dtor.material("Al6061T6");
0035   Material steel_a53  = dtor.material("StainlessSteelA53");
0036
0037   // build magnet components
0038   buildTubeElement(sdet, assembly, dtor, x_det, "yoke", iron);
0039   buildTubeElement(sdet, assembly, dtor, x_det, "coil", nbti);
0040   buildTubeElement(sdet, assembly, dtor, x_det, "tube", steel_304l);
0041   buildPolyElement(sdet, assembly, dtor, x_det, "endplate", steel_304l);
0042   buildTubeElement(sdet, assembly, dtor, x_det, "yokeshield", steel_304l);
0043   buildTubeElement(sdet, assembly, dtor, x_det, "heatshieldbarrel", alum);
0044   buildPolyElement(sdet, assembly, dtor, x_det, "heatshieldend", alum);
0045   buildTubeElement(sdet, assembly, dtor, x_det, "cryobarrel", steel_a53);
0046   buildPolyElement(sdet, assembly, dtor, x_det, "cryoend", steel_a53);
0047
0048   // final placement
0049   auto pv_assembly = dtor.pickMotherVolume(sdet).placeVolume(assembly);
0050   pv_assembly.addPhysVolID("system", det_id);
0051   sdet.setPlacement(pv_assembly);
0052
0053   // update bounding box
0054   assembly->GetShape()->ComputeBBox();
0055
0056   return sdet;
0057 }
0058
0059 void buildPolyElement(dd4hep::DetElement& sdet, dd4hep::Assembly& assembly, dd4hep::Detector& dtor,
0060                       dd4hep::xml::DetElement x_det, string element, dd4hep::Material material) {
0061   // get all elems
0062   xml_coll_t elems_c(x_det, element.c_str());
0063   int elem_id = 1;
0064
0065   // loop over elems
0066   for (; elems_c; ++elems_c) {
0067     // get one element
0068     xml_comp_t elem_c = elems_c;
0069     string elem_name  = elem_c.nameStr();
0070     // get placement coordinates
0071     xml_dim_t elem_pos = elem_c.child(_U(placement));
0072     double elem_theta  = elem_pos.attr<double>("theta");
0073     std::vector<double> z;
0074     std::vector<double> rmax;
0075     std::vector<double> rmin;
0076     // loop over z-planes
0077     xml_coll_t zplanes_c(elem_c, _Unicode(zplane));
0078     for (; zplanes_c; ++zplanes_c) {
0079       // get z-plane
0080       xml_comp_t zplane_c = zplanes_c;
0081       z.push_back(zplane_c.attr<double>(_Unicode(z)));
0082       rmin.push_back(zplane_c.attr<double>(_Unicode(rmin)));
0083       rmax.push_back(zplane_c.attr<double>(_Unicode(rmax)));
0084     }
0085
0086     // set attributes
0087     const string elem_vis =
0088         dd4hep::getAttrOrDefault<std::string>(elem_c, _Unicode(vis), "FFMagnetVis");
0089     sdet.setAttributes(dtor, assembly, x_det.regionStr(), x_det.limitsStr(), elem_vis);
0090
0091     // build solid
0092     Polycone elem_tube(0, 2.0 * M_PI, rmin, rmax, z);
0093     Solid elem_final(elem_tube);
0094
0095     // get all adds
0096     xml_coll_t adds_c(elem_c, _Unicode(add));
0097     // loop over adds
0098     for (; adds_c; ++adds_c) {
0099       Solid add_elem;
0100       // get one cut
0101       xml_comp_t add_c = adds_c;
0102       // get shape
0103       string add_shape = add_c.attr<string>("shape");
0104       // get placement coordinates
0105       xml_dim_t add_pos = add_c.child(_U(placement));
0106       double add_theta  = add_pos.attr<double>("theta");
0107       // get dimentions
0108       xml_dim_t add_dim = add_c.child(_U(dimensions));
0109
0110       Transform3D tf_tmp(RotationZYX(0, add_theta, 0),
0111                          Position(add_pos.x(), add_pos.y(), add_pos.z()));
0112
0113       if (add_shape == "Prism") {
0114         double add_pdx1 = add_dim.attr<double>("pdx1");
0115         double add_pdx2 = add_dim.attr<double>("pdx2");
0116         double add_pdy1 = add_dim.attr<double>("pdy1");
0117         double add_pdy2 = add_dim.attr<double>("pdy2");
0118         double add_pdz  = add_dim.attr<double>("pdz");
0119
0120         // build a solid
0121         Trapezoid add_prism(add_pdx1, add_pdx2, add_pdy1, add_pdy2, add_pdz);
0122         add_elem = add_prism;
0123       } else if (add_shape == "Tube") {
0124         double add_rmin   = add_dim.attr<double>("rmin");
0125         double add_rmax   = add_dim.attr<double>("rmax");
0126         double add_half_l = add_dim.attr<double>("half_length");
0127
0128         // build a solid
0129         Tube add_tube(add_rmin, add_rmax, add_half_l);
0130         add_elem = add_tube;
0131       }
0132
0133       // unite the add with the element solid
0134       elem_final = UnionSolid("elem_final", elem_final, add_elem, tf_tmp);
0135     }
0136
0137     // get all cuts
0138     xml_coll_t cuts_c(elem_c, _Unicode(cut));
0139
0140     Solid cut_final;
0141     // loop over cuts
0142     for (; cuts_c; ++cuts_c) {
0143       // get one cut
0144       xml_comp_t cut_c = cuts_c;
0145       // get shape
0146       string cut_shape = dd4hep::getAttrOrDefault<std::string>(cut_c, _Unicode(shape), "Tube");
0147
0148       // get placement coordinates
0149       xml_dim_t cut_pos = cut_c.child(_U(placement));
0150       double cut_rotX   = cut_pos.attr<double>("rotX");
0151       double cut_rotY   = cut_pos.attr<double>("rotY");
0152       double cut_rotZ   = cut_pos.attr<double>("rotZ");
0153       // get rotation coordinates
0154       xml_dim_t cut_rot    = cut_c.child(_U(rotation));
0155       int cut_rot_num      = cut_rot.attr<int>("num");
0156       double cut_rot_step  = cut_rot.attr<double>("step");
0157       double cut_rot_start = cut_rot.attr<double>("start");
0158       string cut_rot_axis  = cut_rot.attr<string>("axis");
0159
0160       if (cut_shape == "Cone") {
0161         // get dimentions
0162         xml_dim_t cut_dim = cut_c.child(_U(dimensions));
0163         double cut_rmin1  = cut_dim.attr<double>("rmin1");
0164         double cut_rmax1  = cut_dim.attr<double>("rmax1");
0165         double cut_rmin2  = cut_dim.attr<double>("rmin2");
0166         double cut_rmax2  = cut_dim.attr<double>("rmax2");
0167         double cut_dz     = cut_dim.attr<double>("dz");
0168
0169         // build a solid
0170         Cone cut_cone(cut_dz, cut_rmin1, cut_rmax1, cut_rmin2, cut_rmax2);
0171         cut_final = cut_cone;
0172       } else if (cut_shape == "Tube") {
0173         // get dimentions
0174         xml_dim_t cut_dim = cut_c.child(_U(dimensions));
0175         double cut_rmin   = cut_dim.attr<double>("rmin");
0176         double cut_rmax   = cut_dim.attr<double>("rmax");
0177         double cut_half_l = cut_dim.attr<double>("half_length");
0178
0179         // build a solid
0180         Tube cut_tube(cut_rmin, cut_rmax, cut_half_l);
0181         cut_final = cut_tube;
0182       }
0183
0184       // loop over rot steps
0185       for (int i = 0; i < cut_rot_num; i++) {
0186         Position pos_tmp(cut_pos.x(), cut_pos.y(), cut_pos.z());
0187         double ang_tmp = cut_rot_start + i * cut_rot_step;
0188         Rotation3D rot_tmp;
0189         if (cut_rot_axis == "X") {
0190           rot_tmp = RotationX(ang_tmp);
0191         } else if (cut_rot_axis == "Y") {
0192           rot_tmp = RotationY(ang_tmp);
0193         } else {
0194           rot_tmp = RotationZ(ang_tmp);
0195         }
0196         pos_tmp = rot_tmp * pos_tmp;
0197
0198         Transform3D tf_tmp(RotationZYX(cut_rotZ, cut_rotY, cut_rotX), pos_tmp);
0199         // subtract the cut from the element solid
0200         elem_final = SubtractionSolid("elem_final", elem_final, cut_final, tf_tmp);
0201       }
0202     }
0203
0204     // create volume
0205     Volume elem_vol(elem_name, elem_final, material);
0206
0207     // placement
0208     auto elem_pv = assembly.placeVolume(
0209         elem_vol, Transform3D(Translation3D(elem_pos.x(), elem_pos.y(), elem_pos.z()) *
0210                               RotationY(elem_theta)));
0211     elem_pv.addPhysVolID(element, elem_id);
0212     DetElement elem_de(sdet, elem_name, elem_id);
0213     elem_de.setPlacement(elem_pv);
0214     elem_de.setAttributes(dtor, elem_vol, x_det.regionStr(), x_det.limitsStr(), elem_vis);
0215     elem_id++;
0216   }
0217
0218   return;
0219 }
0220
0221 void buildTubeElement(dd4hep::DetElement& sdet, dd4hep::Assembly& assembly, dd4hep::Detector& dtor,
0222                       dd4hep::xml::DetElement x_det, string element, dd4hep::Material material) {
0223   // get all elems
0224   xml_coll_t elems_c(x_det, element.c_str());
0225   int elem_id = 1;
0226
0227   // loop over elems
0228   for (; elems_c; ++elems_c) {
0229     // get one element
0230     xml_comp_t elem_c = elems_c;
0231     string elem_name  = elem_c.nameStr();
0232     // get placement coordinates
0233     xml_dim_t elem_pos = elem_c.child(_U(placement));
0234     double elem_theta  = elem_pos.attr<double>("theta");
0235
0236     // set attributes
0237     const string elem_vis =
0238         dd4hep::getAttrOrDefault<std::string>(elem_c, _Unicode(vis), "FFMagnetVis");
0239     sdet.setAttributes(dtor, assembly, x_det.regionStr(), x_det.limitsStr(), elem_vis);
0240
0241     // get shape
0242     string elem_shape = dd4hep::getAttrOrDefault<std::string>(elem_c, _Unicode(shape), "Tube");
0243
0244     Solid elem_sub;
0245     if (elem_shape == "Cone") {
0246       // get dimentions
0247       xml_dim_t elem_dim = elem_c.child(_U(dimensions));
0248       double elem_rmin1  = elem_dim.attr<double>("rmin1");
0249       double elem_rmax1  = elem_dim.attr<double>("rmax1");
0250       double elem_rmin2  = elem_dim.attr<double>("rmin2");
0251       double elem_rmax2  = elem_dim.attr<double>("rmax2");
0252       double elem_dz     = elem_dim.attr<double>("dz");
0253
0254       // build a solid
0255       Cone elem_cone(elem_dz, elem_rmin1, elem_rmax1, elem_rmin2, elem_rmax2);
0256       elem_sub = elem_cone;
0257     } else if (elem_shape == "ConeSegment") {
0258       // get dimentions
0259       xml_dim_t elem_dim = elem_c.child(_U(dimensions));
0260       double elem_rmin1  = elem_dim.attr<double>("rmin1");
0261       double elem_rmax1  = elem_dim.attr<double>("rmax1");
0262       double elem_rmin2  = elem_dim.attr<double>("rmin2");
0263       double elem_rmax2  = elem_dim.attr<double>("rmax2");
0264       double elem_dz     = elem_dim.attr<double>("dz");
0265       double elem_sphi   = elem_dim.attr<double>("sphi");
0266       double elem_dphi   = elem_dim.attr<double>("dphi");
0267
0268       // build a solid
0269       ConeSegment elem_conesegment(elem_dz, elem_rmin1, elem_rmax1, elem_rmin2, elem_rmax2,
0270                                    elem_sphi, elem_sphi + elem_dphi);
0271       elem_sub = elem_conesegment;
0272     } else if (elem_shape == "Tube") {
0273       // get dimentions
0274       xml_dim_t elem_dim = elem_c.child(_U(dimensions));
0275       double elem_rmin   = elem_dim.attr<double>("rmin");
0276       double elem_rmax   = elem_dim.attr<double>("rmax");
0277       double elem_half_l = elem_dim.attr<double>("half_length");
0278       double elem_sphi   = elem_dim.attr<double>("sphi");
0279       double elem_dphi   = elem_dim.attr<double>("dphi");
0280
0281       // build solid
0282       Tube elem_tube(elem_rmin, elem_rmax, elem_half_l, elem_sphi, elem_sphi + elem_dphi);
0283       elem_sub = elem_tube;
0284     }
0285
0286     Solid elem_final(elem_sub);
0287
0288     // combine sub-elements
0289     if (elem_pos.hasAttr(_Unicode(phiNum))) {
0290       int phi_num      = elem_pos.attr<int>("phiNum");
0291       double phi_step  = elem_pos.attr<double>("phiStep");
0292       double phi_start = elem_pos.attr<double>("phiStart");
0293
0294       // loop over steps
0295       for (int i = 0; i < phi_num; i++) {
0296         double phi_tmp = phi_start + i * phi_step;
0297         Transform3D tf_tmp(RotationZ(phi_tmp), Position(0, 0, 0));
0298         // unite sub-elements
0299         elem_final = UnionSolid("elem_final", elem_final, elem_sub, tf_tmp);
0300       }
0301     }
0302
0303     // get all cuts
0304     xml_coll_t cuts_c(elem_c, _Unicode(cut));
0305     // loop over cuts
0306     for (; cuts_c; ++cuts_c) {
0307       // get one cut
0308       xml_comp_t cut_c = cuts_c;
0309       // get shape
0310       string add_shape = dd4hep::getAttrOrDefault<std::string>(cut_c, _Unicode(shape), "Tube");
0311       // get placement coordinates
0312       xml_dim_t cut_pos = cut_c.child(_U(placement));
0313       double cut_rotX   = cut_pos.attr<double>("rotX");
0314       double cut_rotY   = cut_pos.attr<double>("rotY");
0315       double cut_rotZ   = cut_pos.attr<double>("rotZ");
0316       // get rotation coordinates
0317       xml_dim_t cut_rot    = cut_c.child(_U(rotation));
0318       int cut_rot_num      = cut_rot.attr<int>("num");
0319       double cut_rot_step  = cut_rot.attr<double>("step");
0320       double cut_rot_start = cut_rot.attr<double>("start");
0321       string cut_rot_axis  = cut_rot.attr<string>("axis");
0322
0323       Solid cut_elem;
0324       if (add_shape == "Box") {
0325         // get dimentions
0326         xml_dim_t cut_dim = cut_c.child(_U(dimensions));
0327         double cut_dx     = cut_dim.attr<double>("dx");
0328         double cut_dy     = cut_dim.attr<double>("dy");
0329         double cut_dz     = cut_dim.attr<double>("dz");
0330
0331         // build a solid
0332         Box cut_box(cut_dx, cut_dy, cut_dz);
0333         cut_elem = cut_box;
0334       } else if (add_shape == "Tube") {
0335         // get dimentions
0336         xml_dim_t cut_dim = cut_c.child(_U(dimensions));
0337         double cut_rmin   = cut_dim.attr<double>("rmin");
0338         double cut_rmax   = cut_dim.attr<double>("rmax");
0339         double cut_half_l = cut_dim.attr<double>("half_length");
0340
0341         // build a solid
0342         Tube cut_tube(cut_rmin, cut_rmax, cut_half_l);
0343         cut_elem = cut_tube;
0344       } else if (add_shape == "Cone") {
0345         // get dimentions
0346         xml_dim_t cut_dim = cut_c.child(_U(dimensions));
0347         double cut_rmin1  = cut_dim.attr<double>("rmin1");
0348         double cut_rmax1  = cut_dim.attr<double>("rmax1");
0349         double cut_rmin2  = cut_dim.attr<double>("rmin2");
0350         double cut_rmax2  = cut_dim.attr<double>("rmax2");
0351         double cut_dz     = cut_dim.attr<double>("dz");
0352
0353         // build a solid
0354         Cone cut_cone(cut_dz, cut_rmin1, cut_rmax1, cut_rmin2, cut_rmax2);
0355         cut_elem = cut_cone;
0356       }
0357       // loop over rot steps
0358       for (int i = 0; i < cut_rot_num; i++) {
0359         Position pos_tmp(cut_pos.x(), cut_pos.y(), cut_pos.z());
0360         double ang_tmp = cut_rot_start + i * cut_rot_step;
0361         Rotation3D rot_tmp;
0362         if (cut_rot_axis == "X") {
0363           rot_tmp = RotationX(ang_tmp);
0364         } else if (cut_rot_axis == "Y") {
0365           rot_tmp = RotationY(ang_tmp);
0366         } else {
0367           rot_tmp = RotationZ(ang_tmp);
0368         }
0369         pos_tmp = rot_tmp * pos_tmp;
0370
0371         Transform3D tf_tmp(RotationZYX(cut_rotZ, cut_rotY, cut_rotX), pos_tmp);
0372         // subtract the cut from the element solid
0373         elem_final = SubtractionSolid("elem_final", elem_final, cut_elem, tf_tmp);
0374       }
0375     }
0376
0377     // create volume
0378     Volume elem_vol(elem_name, elem_final, material);
0379
0380     // placement
0381     auto elem_pv = assembly.placeVolume(
0382         elem_vol, Transform3D(Translation3D(elem_pos.x(), elem_pos.y(), elem_pos.z()) *
0383                               RotationY(elem_theta)));
0384     elem_pv.addPhysVolID(element, elem_id);
0385     DetElement elem_de(sdet, elem_name, elem_id);
0386     elem_de.setPlacement(elem_pv);
0387     elem_de.setAttributes(dtor, elem_vol, x_det.regionStr(), x_det.limitsStr(), elem_vis);
0388     elem_id++;
0389   }
0390
0391   return;
0392 }
0393
0394 DECLARE_DETELEMENT(ip6_CryostatMagnet, build_magnet)