epic/src/ScFiCalorimeter_geo.cpp

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