athena/src/ShashlikCalorimeter_geo.cpp

0001 //==========================================================================
0002 //  Implementation for shashlik calorimeter modules
0003 //  it supports disk placements with (rmin, rmax), and (phimin, phimax)
0004 //--------------------------------------------------------------------------
0005 //  Author: Chao Peng (ANL)
0006 //  Date: 06/22/2021
0007 //==========================================================================
0008
0009 #include "GeometryHelpers.h"
0010 #include "DD4hep/DetFactoryHelper.h"
0011 #include <XML/Layering.h>
0012 #include <XML/Helper.h>
0013 #include <iostream>
0014 #include <algorithm>
0015 #include <tuple>
0016 #include <math.h>
0017
0018 using namespace dd4hep;
0019
0020 static void add_disk_shashlik(Detector& desc, Assembly &env, xml::Collection_t &plm, SensitiveDetector &sens, int id);
0021
0022 // helper function to get x, y, z if defined in a xml component
0023 template<class XmlComp>
0024 Position get_xml_xyz(XmlComp &comp, dd4hep::xml::Strng_t name)
0025 {
0026     Position pos(0., 0., 0.);
0027     if (comp.hasChild(name)) {
0028         auto child = comp.child(name);
0029         pos.SetX(dd4hep::getAttrOrDefault<double>(child, _Unicode(x), 0.));
0030         pos.SetY(dd4hep::getAttrOrDefault<double>(child, _Unicode(y), 0.));
0031         pos.SetZ(dd4hep::getAttrOrDefault<double>(child, _Unicode(z), 0.));
0032     }
0033     return pos;
0034 }
0035
0036 static Ref_t create_detector(Detector& desc, xml::Handle_t handle, SensitiveDetector sens)
0037 {
0038     static const std::string func = "ShashlikCalorimeter";
0039     xml::DetElement detElem = handle;
0040     std::string detName = detElem.nameStr();
0041     int detID = detElem.id();
0042     DetElement det(detName, detID);
0043     sens.setType("calorimeter");
0044     // envelope
0045     Assembly assembly(detName);
0046
0047     // module placement
0048     xml::Component plm = detElem.child(_Unicode(placements));
0049     int sector = 1;
0050     for (xml::Collection_t mod(plm, _Unicode(disk)); mod; ++mod) {
0051         add_disk_shashlik(desc, assembly, mod, sens, sector++);
0052     }
0053
0054     // detector position and rotation
0055     auto pos = get_xml_xyz(detElem, _Unicode(position));
0056     auto rot = get_xml_xyz(detElem, _Unicode(rotation));
0057     Volume motherVol = desc.pickMotherVolume(det);
0058     Transform3D tr = Translation3D(pos.x(), pos.y(), pos.z()) * RotationZYX(rot.z(), rot.y(), rot.x());
0059     PlacedVolume envPV = motherVol.placeVolume(assembly, tr);
0060     envPV.addPhysVolID("system", detID);
0061     det.setPlacement(envPV);
0062     return det;
0063 }
0064
0065 // helper function to build module with or w/o wrapper
0066 std::tuple<Volume, int, double, double> build_shashlik(Detector &desc, xml::Collection_t &plm, SensitiveDetector &sens)
0067 {
0068     auto mod = plm.child(_Unicode(module));
0069     // a modular volume
0070     std::string shape = dd4hep::getAttrOrDefault(mod, _Unicode(shape), "square");
0071     std::transform(shape.begin(), shape.end(), shape.begin(), [](char c) { return std::tolower(c); });
0072     int nsides = 4;
0073     if (shape == "hexagon") {
0074         nsides = 6;
0075     } else if (shape != "square") {
0076         std::cerr << "ShashlikCalorimeter Error: Unsupported shape of module " << shape
0077                   << ". Please choose from (square, hexagon). Proceed with square shape." << std::endl;
0078     }
0079     double slen = mod.attr<double>(_Unicode(side_length));
0080     double rmax = slen/2./std::sin(M_PI/nsides);
0081     Layering layering(mod);
0082     auto len = layering.totalThickness();
0083
0084     // wrapper info
0085     PolyhedraRegular mpoly(nsides, 0., rmax, len);
0086     Volume mvol("shashlik_module_vol", mpoly, desc.air());
0087     mvol.setVisAttributes(desc.visAttributes(dd4hep::getAttrOrDefault(mod, _Unicode(vis), "GreenVis")));
0088
0089     double gap = 0.;
0090     Volume wvol("shashlik_wrapper_vol");
0091     if (plm.hasChild(_Unicode(wrapper))) {
0092         auto wrap = plm.child(_Unicode(wrapper));
0093         gap = wrap.attr<double>(_Unicode(thickness));
0094         if (gap > 1e-6*mm) {
0095             wvol.setSolid(PolyhedraRegular(nsides, 0., rmax + gap, len));
0096             wvol.setMaterial(desc.material(wrap.attr<std::string>(_Unicode(material))));
0097             wvol.setVisAttributes(desc.visAttributes(dd4hep::getAttrOrDefault(wrap, _Unicode(vis), "WhiteVis")));
0098             wvol.placeVolume(mvol, Position{0., 0., 0.});
0099         }
0100     }
0101
0102     // layer start point
0103     double lz  = -len/2.;
0104     int lnum = 1;
0105     // Loop over the sets of layer elements in the detector.
0106     for (xml_coll_t li(mod, _U(layer)); li; ++li) {
0107         int repeat = li.attr<int>(_Unicode(repeat));
0108         // Loop over number of repeats for this layer.
0109         for (int j = 0; j < repeat; j++) {
0110             std::string lname = Form("layer%d", lnum);
0111             double lthick = layering.layer(lnum - 1)->thickness();  // Layer's thickness.
0112             PolyhedraRegular lpoly(nsides, 0., rmax, lthick);
0113             Volume lvol(lname, lpoly, desc.air());
0114
0115             // Loop over the sublayers or slices for this layer.
0116             int snum = 1;
0117             double sz = -lthick/2.;
0118             for (xml_coll_t si(li, _U(slice)); si; ++si) {
0119                 std::string sname = Form("slice%d", snum);
0120                 double sthick = si.attr<double>(_Unicode(thickness));
0121                 PolyhedraRegular spoly(nsides, 0., rmax, sthick);
0122                 Volume svol(sname, spoly, desc.material(si.attr<std::string>(_Unicode(material))));
0123
0124                 std::string issens = dd4hep::getAttrOrDefault(si, _Unicode(sensitive), "no");
0125                 std::transform(issens.begin(), issens.end(), issens.begin(), [](char c) { return std::tolower(c); });
0126                 if ((issens == "yes") || (issens == "y") || (issens == "true")) {
0127                     svol.setSensitiveDetector(sens);
0128                 }
0129                 svol.setAttributes(desc,
0130                         dd4hep::getAttrOrDefault(si, _Unicode(region), ""),
0131                         dd4hep::getAttrOrDefault(si, _Unicode(limits), ""),
0132                         dd4hep::getAttrOrDefault(si, _Unicode(vis), "InvisibleNoDaughters"));
0133
0134                 // Slice placement.
0135                 auto slicePV = lvol.placeVolume(svol, Position(0, 0, sz + sthick/2.));
0136                 slicePV.addPhysVolID("slice", snum++);
0137                 // Increment Z position of slice.
0138                 sz += sthick;
0139             }
0140
0141             // Set region, limitset, and vis of layer.
0142             lvol.setAttributes(desc,
0143                     dd4hep::getAttrOrDefault(li, _Unicode(region), ""),
0144                     dd4hep::getAttrOrDefault(li, _Unicode(limits), ""),
0145                     dd4hep::getAttrOrDefault(li, _Unicode(vis), "InvisibleNoDaughters"));
0146
0147             auto layerPV = mvol.placeVolume(lvol, Position(0, 0, lz + lthick/2));
0148             layerPV.addPhysVolID("layer", lnum++);
0149             // Increment to next layer Z position.
0150             lz += lthick;
0151         }
0152     }
0153
0154     if (gap > 1e-6*mm) {
0155         return std::make_tuple(wvol, nsides, 2.*std::sin(M_PI/nsides)*(rmax + gap), len);
0156     } else {
0157         return std::make_tuple(mvol, nsides, slen, len);
0158     }
0159 }
0160
0161 // place disk of modules
0162 static void add_disk_shashlik(Detector& desc, Assembly &env, xml::Collection_t &plm, SensitiveDetector &sens, int sid)
0163 {
0164     auto [mvol, nsides, sidelen, len]  = build_shashlik(desc, plm, sens);
0165     int sector_id = dd4hep::getAttrOrDefault<int>(plm, _Unicode(sector), sid);
0166     int id_begin = dd4hep::getAttrOrDefault<int>(plm, _Unicode(id_begin), 1);
0167     double rmin = plm.attr<double>(_Unicode(rmin));
0168     double rmax = plm.attr<double>(_Unicode(rmax));
0169     double phimin = dd4hep::getAttrOrDefault<double>(plm, _Unicode(phimin), 0.);
0170     double phimax = dd4hep::getAttrOrDefault<double>(plm, _Unicode(phimax), 2.*M_PI);
0171
0172     auto points = (nsides == 6) ? athena::geo::fillHexagons({0., 0.}, sidelen, rmin, rmax, phimin, phimax)
0173                                 : athena::geo::fillSquares({0., 0.}, sidelen*1.414, rmin, rmax, phimin, phimax);
0174     // placement to mother
0175     auto pos = get_xml_xyz(plm, _Unicode(position));
0176     auto rot = get_xml_xyz(plm, _Unicode(rotation));
0177     int mid = 0;
0178     for (auto &p : points) {
0179         Transform3D tr = RotationZYX(rot.z(), rot.y(), rot.x())
0180                        * Translation3D(pos.x() + p.x(), pos.y() + p.y(), pos.z() + len/2.)
0181                        * RotationZ((nsides == 4) ? 45*degree : 0);
0182         auto modPV = env.placeVolume(mvol, tr);
0183         modPV.addPhysVolID("sector", sector_id).addPhysVolID("module", id_begin + mid++);
0184     }
0185 }
0186
0187
0188 DECLARE_DETELEMENT(ShashlikCalorimeter, create_detector)
0189