athena/src/ScFiCalorimeter_geo.cpp

0001 //==========================================================================
0002 //  Scintillating fiber calorimeter with tower shape blocks
0003 //  reference: https://github.com/adamjaro/lmon/blob/master/calo/src/WScFiZXv3.cxx
0004 //  Support disk placement
0005 //--------------------------------------------------------------------------
0006 //  Author: Chao Peng (ANL)
0007 //  Date: 07/19/2021
0008 //==========================================================================
0009
0010 #include "GeometryHelpers.h"
0011 #include "DD4hep/DetFactoryHelper.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 using Point =  ROOT::Math::XYPoint;
0020
0021 std::tuple<Volume, Position> build_module(const Detector &desc, const xml::Component &mod_x, SensitiveDetector &sens);
0022
0023 // helper function to get x, y, z if defined in a xml component
0024 template<class XmlComp>
0025 Position get_xml_xyz(const XmlComp &comp, dd4hep::xml::Strng_t name)
0026 {
0027     Position pos(0., 0., 0.);
0028     if (comp.hasChild(name)) {
0029         auto child = comp.child(name);
0030         pos.SetX(dd4hep::getAttrOrDefault<double>(child, _Unicode(x), 0.));
0031         pos.SetY(dd4hep::getAttrOrDefault<double>(child, _Unicode(y), 0.));
0032         pos.SetZ(dd4hep::getAttrOrDefault<double>(child, _Unicode(z), 0.));
0033     }
0034     return pos;
0035 }
0036
0037 // main
0038 static Ref_t create_detector(Detector& desc, xml::Handle_t handle, SensitiveDetector sens)
0039 {
0040     xml::DetElement detElem = handle;
0041     std::string detName = detElem.nameStr();
0042     int detID = detElem.id();
0043     DetElement det(detName, detID);
0044     sens.setType("calorimeter");
0045     auto dim = detElem.dimensions();
0046     auto rmin = dim.rmin();
0047     auto rmax = dim.rmax();
0048     auto length = dim.length();
0049     auto phimin = dd4hep::getAttrOrDefault<double>(dim, _Unicode(phimin), 0.);
0050     auto phimax = dd4hep::getAttrOrDefault<double>(dim, _Unicode(phimax), 2.*M_PI);
0051     // envelope
0052     Tube envShape(rmin, rmax, length/2., phimin, phimax);
0053     Volume env(detName + "_envelope", envShape, desc.material("Air"));
0054     env.setVisAttributes(desc.visAttributes(detElem.visStr()));
0055
0056     // build module
0057     auto [modVol, modSize] = build_module(desc, detElem.child(_Unicode(module)), sens);
0058     double modSizeR = std::sqrt(modSize.x() * modSize.x() + modSize.y() * modSize.y());
0059     double assembly_rwidth = modSizeR*2.;
0060     int nas = int((rmax - rmin) / assembly_rwidth) + 1;
0061     std::vector<Assembly> assemblies;
0062     // calorimeter block z-offsets (as blocks are shorter than the volume length)
0063     const double block_offset = -0.5*(length - modSize.z());
0064     for (int i = 0; i < nas; ++i) {
0065         Assembly assembly(detName + Form("_ring%d", i + 1));
0066         auto assemblyPV = env.placeVolume(assembly, Position{0., 0., block_offset});
0067         assemblyPV.addPhysVolID("ring", i + 1);
0068         assemblies.emplace_back(std::move(assembly));
0069     }
0070     // std::cout << assemblies.size() << std::endl;
0071
0072     int modid = 1;
0073     for (int ix = 0; ix < int(2.*rmax / modSize.x()) + 1; ++ix) {
0074         double mx = modSize.x() * ix - rmax;
0075         for (int iy = 0; iy < int(2.*rmax / modSize.y()) + 1; ++iy) {
0076             double my = modSize.y() * iy - rmax;
0077             double mr = std::sqrt(mx*mx + my*my);
0078             if (mr - modSizeR >= rmin && mr + modSizeR <= rmax) {
0079                 int ias = int((mr - rmin) / assembly_rwidth);
0080                 auto &assembly = assemblies[ias];
0081                 auto modPV = assembly.placeVolume(modVol, Position(mx, my, 0.));
0082                 modPV.addPhysVolID("module", modid++);
0083             }
0084         }
0085     }
0086
0087     desc.add(Constant(detName + "_NModules", std::to_string(modid - 1)));
0088
0089     for (auto &assembly : assemblies) {
0090         assembly.ptr()->Voxelize("");
0091     }
0092
0093     // detector position and rotation
0094     auto pos = get_xml_xyz(detElem, _Unicode(position));
0095     auto rot = get_xml_xyz(detElem, _Unicode(rotation));
0096     Volume motherVol = desc.pickMotherVolume(det);
0097     Transform3D tr = Translation3D(pos.x(), pos.y(), pos.z()) * RotationZYX(rot.z(), rot.y(), rot.x());
0098     PlacedVolume envPV = motherVol.placeVolume(env, tr);
0099     envPV.addPhysVolID("system", detID);
0100     det.setPlacement(envPV);
0101     return det;
0102 }
0103
0104 // helper function to build module with scintillating fibers
0105 std::tuple<Volume, Position> build_module(const Detector &desc, const xml::Component &mod_x, SensitiveDetector &sens)
0106 {
0107     auto sx = mod_x.attr<double>(_Unicode(sizex));
0108     auto sy = mod_x.attr<double>(_Unicode(sizey));
0109     auto sz = mod_x.attr<double>(_Unicode(sizez));
0110
0111     Box modShape(sx/2., sy/2., sz/2.);
0112     auto modMat = desc.material(mod_x.attr<std::string>(_Unicode(material)));
0113     Volume modVol("module_vol", modShape, modMat);
0114     if (mod_x.hasAttr(_Unicode(vis))) {
0115         modVol.setVisAttributes(desc.visAttributes(mod_x.attr<std::string>(_Unicode(vis))));
0116     }
0117
0118     if (mod_x.hasChild("fiber")) {
0119       auto fiber_x  = mod_x.child(_Unicode(fiber));
0120       auto fr       = fiber_x.attr<double>(_Unicode(radius));
0121       auto fsx      = fiber_x.attr<double>(_Unicode(spacex));
0122       auto fsy      = fiber_x.attr<double>(_Unicode(spacey));
0123       auto foff     = dd4hep::getAttrOrDefault<double>(fiber_x, _Unicode(offset), 0.5*mm);
0124       auto fiberMat = desc.material(fiber_x.attr<std::string>(_Unicode(material)));
0125       Tube fiberShape(0., fr, sz/2.);
0126       Volume fiberVol("fiber_vol", fiberShape, fiberMat);
0127       fiberVol.setSensitiveDetector(sens);
0128
0129       // Fibers are placed in a honeycomb with the radius = sqrt(3)/2. * hexagon side length
0130       // So each fiber is fully contained in a regular hexagon, which are placed as
0131       //           ______________________________________
0132       //           |          ____        ____          |
0133       // even:     |         /    \      /    \         |
0134       //           |    ____/      \____/      \____    |
0135       //           |   /    \      /    \      /    \   |
0136       // odd:      |  /      \____/      \____/      \  |
0137       //           |  \      /    \      /    \      /  |
0138       //           |   \____/      \____/      \____/   |
0139       // even:     |        \      /    \      /        |
0140       //           |         \____/      \____/      ___|___
0141       //           |____________________________________|___offset
0142       //                                              | |
0143       //                                              |offset
0144       // the parameters space x and space y are used to add additional spaces between the hexagons
0145       double fside  = 2. / std::sqrt(3.) * fr;
0146       double fdistx = 2. * fside + fsx;
0147       double fdisty = 2. * fr + fsy;
0148
0149       // maximum numbers of the fibers, help narrow the loop range
0150       int nx = int(sx / (2.*fr)) + 1;
0151       int ny = int(sy / (2.*fr)) + 1;
0152
0153       // std::cout << sx << ", " << sy << ", " << fr << ", " << nx << ", " << ny << std::endl;
0154
0155       // place the fibers
0156       double y0 = (foff + fside);
0157       int nfibers = 0;
0158       for (int iy = 0; iy < ny; ++iy) {
0159           double y = y0 + fdisty * iy;
0160           // about to touch the boundary
0161           if ((sy - y) < y0) { break; }
0162           double x0 = (iy % 2) ? (foff + fside) : (foff + fside + fdistx / 2.);
0163           for (int ix = 0; ix < nx; ++ix) {
0164               double x = x0 + fdistx * ix;
0165               // about to touch the boundary
0166               if ((sx - x) < x0) { break; }
0167               auto fiberPV = modVol.placeVolume(fiberVol, nfibers++, Position{x - sx/2., y - sy/2., 0});
0168               //std::cout << "(" << ix << ", " << iy << ", " << x - sx/2. << ", " << y - sy/2. << ", " << fr << "),\n";
0169               fiberPV.addPhysVolID("fiber_x", ix + 1).addPhysVolID("fiber_y", iy + 1);
0170           }
0171       }
0172     // if no fibers we make the module itself sensitive
0173     } else {
0174       modVol.setSensitiveDetector(sens);
0175     }
0176
0177
0178     return std::make_tuple(modVol, Position{sx, sy, sz});
0179 }
0180
0181 DECLARE_DETELEMENT(ScFiCalorimeter, create_detector)
0182