Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2024-11-15 08:59:19

0001 //==========================================================================
0002 //  A general implementation for homogeneous calorimeter
0003 //  it supports three types of placements
0004 //  1. Individual module placement with module dimensions and positions
0005 //  2. Array placement with module dimensions and numbers of row and columns
0006 //  3. Disk placement with module dimensions and (Rmin, Rmax), and (Phimin, Phimax)
0007 //  4. Lines placement with module dimensions and (mirrorx, mirrory)
0008 //     (NOTE: anchor point is the 0th block of the line instead of line center)
0009 //--------------------------------------------------------------------------
0010 //  Author: Chao Peng (ANL)
0011 //  Date: 06/09/2021
0012 //==========================================================================
0013 
0014 #include "GeometryHelpers.h"
0015 #include "DD4hep/DetFactoryHelper.h"
0016 #include <XML/Helper.h>
0017 #include <iostream>
0018 #include <algorithm>
0019 #include <tuple>
0020 #include <math.h>
0021 
0022 using namespace dd4hep;
0023 
0024 /** \addtogroup calorimeters Calorimeters
0025  */
0026 
0027 /** \addtogroup Homogeneous Calorimeter
0028  * \brief Type: **HomogeneousCalorimeter**.
0029  * \author C. Peng
0030  * \ingroup calorimeters
0031  *
0032  *
0033  * \code
0034  *   <detector id="1" name="HyCal" type="HomogeneousCalorimeter" readout="EcalHits">
0035  *     <position x="0" y="0" z="0"/>
0036  *     <rotation x="0" y="0" z="0"/>
0037  *     <placements>
0038  *       <array nrow="34" ncol="34" sector="1">
0039  *         <position x="0" y="0" z="-9.73*cm"/>
0040  *         <module sizex="2.05*cm" sizey="2.05*cm" sizez="18*cm" vis="GreenVis" material="PbWO4"/>
0041  *         <wrapper thickness="0.015*cm" material="Epoxy" vis="WhiteVis"/>
0042  *         <removal row="16" col="16"/>
0043  *         <removal row="16" col="17"/>
0044  *         <removal row="17" col="16"/>
0045  *         <removal row="17" col="17"/>
0046  *       </array>
0047  *       <array nrow="6" ncol="24" sector="2">
0048  *         <position x="-17*(2.05+0.015)*cm+12*(3.8+0.015)*cm" y="17*(2.05+0.015)*cm+3*(3.8+0.015)*cm" z="0"/>
0049  *         <module sizex="3.8*cm" sizey="3.8*cm" sizez="45*cm" vis="BlueVis" material="PbGlass"/>
0050  *         <wrapper thickness="0.015*cm" material="Epoxy" vis="WhiteVis"/>
0051  *       </array>
0052  *       <array nrow="24" ncol="6" sector="3">
0053  *         <position x="17*(2.05+0.015)*cm+3*(3.8+0.015)*cm" y="17*(2.05+0.015)*cm-12*(3.8+0.015)*cm" z="0"/>
0054  *         <module sizex="3.8*cm" sizey="3.8*cm" sizez="45*cm" vis="BlueVis" material="PbGlass"/>
0055  *         <wrapper thickness="0.015*cm" material="Epoxy" vis="WhiteVis"/>
0056  *       </array>
0057  *       <array nrow="6" ncol="24" sector="4">
0058  *         <position x="17*(2.05+0.015)*cm-12*(3.8+0.015)*cm" y="-17*(2.05+0.015)*cm-3*(3.8+0.015)*cm" z="0"/>
0059  *         <module sizex="3.8*cm" sizey="3.8*cm" sizez="45*cm" vis="BlueVis" material="PbGlass"/>
0060  *         <wrapper thickness="0.015*cm" material="Epoxy" vis="WhiteVis"/>
0061  *       </array>
0062  *       <array nrow="24" ncol="6" sector="3">
0063  *         <position x="-17*(2.05+0.015)*cm-3*(3.8+0.015)*cm" y="-17*(2.05+0.015)*cm+12*(3.8+0.015)*cm" z="0"/>
0064  *         <module sizex="3.8*cm" sizey="3.8*cm" sizez="45*cm" vis="BlueVis" material="PbGlass"/>
0065  *         <wrapper thickness="0.015*cm" material="Epoxy" vis="WhiteVis"/>
0066  *       </array>
0067  *     </placements>
0068  *   </detector>
0069  *
0070  *   <detector id="2" name="SomeBlocks" type="HomogeneousCalorimeter" readout="EcalHits">
0071  *     <position x="0" y="0" z="30*cm"/>
0072  *     <rotation x="0" y="0" z="0"/>
0073  *     <placements>
0074  *       <individuals sector="1"/>
0075  *         <module sizex="2.05*cm" sizey="2.05*cm" sizez="20*cm" vis="GreenVis" material="PbWO4"/>
0076  *         <wrapper thickness="0.015*cm" material="Epoxy" vis="WhiteVis"/>
0077  *         <placement x="1*cm" y="1*cm" z="0" id="1"/>
0078  *         <placement x="-1*cm" y="1*cm" z="0" id="2"/>
0079  *         <placement x="1*cm" y="-1*cm" z="0" id="3"/>
0080  *         <placement x="-1*cm" y="-1*cm" z="0" id="4"/>
0081  *       </individuals>
0082  *     </placements>
0083  *   </detector>
0084  *
0085  *   <detector id="2" name="DiskShapeCalorimeter" type="HomogeneousCalorimeter" readout="EcalHits">
0086  *     <position x="0" y="0" z="-30*cm"/>
0087  *     <rotation x="0" y="0" z="0"/>
0088  *     <placements>
0089  *       <disk rmin="25*cm" rmax="125*cm" length="20.5*cm" phimin="0" phimax="360*degree" sector="1"/>
0090  *         <module sizex="2.05*cm" sizey="2.05*cm" sizez="20*cm" vis="GreenVis" material="PbWO4"/>
0091  *         <wrapper thickness="0.015*cm" material="Epoxy" vis="WhiteVis"/>
0092  *     </placements>
0093  *   </detector>
0094  *
0095  *   <detector id="3" name="SomeLines" type="HomogeneousCalorimeter" readout="EcalHits">
0096  *     <position x="0" y="0" z="60*cm"/>
0097  *     <rotation x="0" y="0" z="0"/>
0098  *     <placements>
0099  *       <lines sector="1" mirrorx="true" mirrory="true"/>
0100  *         <module sizex="2.05*cm" sizey="2.05*cm" sizez="20*cm" vis="GreenVis" material="PbWO4"/>
0101  *         <wrapper thickness="0.015*cm" material="Epoxy" vis="WhiteVis"/>
0102  *         <line x="10.25*mm" y="10.25*mm" axis="x" begin="8" nmods="16"/>
0103  *         <line x="10.25*mm" y="30.75*mm" axis="y" begin="8" nmods="16"/>
0104  *         <line x="10.25*mm" y="51.25*mm" axis="z" begin="8" nmods="16"/>
0105  *       </individuals>
0106  *     </placements>
0107  *   </detector>
0108  * \endcode
0109  *
0110  * @{
0111  */
0112 
0113 // headers
0114 static std::tuple<int, int> add_individuals(Detector& desc, Assembly &env, xml::Collection_t &plm, SensitiveDetector &sens, int id);
0115 static std::tuple<int, int> add_array(Detector& desc, Assembly &env, xml::Collection_t &plm, SensitiveDetector &sens, int id);
0116 static std::tuple<int, int> add_disk(Detector& desc, Assembly &env, xml::Collection_t &plm, SensitiveDetector &sens, int id);
0117 static std::tuple<int, int> add_lines(Detector& desc, Assembly &env, xml::Collection_t &plm, SensitiveDetector &sens, int id);
0118 
0119 // helper function to get x, y, z if defined in a xml component
0120 template<class XmlComp>
0121 Position get_xml_xyz(XmlComp &comp, dd4hep::xml::Strng_t name)
0122 {
0123     Position pos(0., 0., 0.);
0124     if (comp.hasChild(name)) {
0125         auto child = comp.child(name);
0126         pos.SetX(dd4hep::getAttrOrDefault<double>(child, _Unicode(x), 0.));
0127         pos.SetY(dd4hep::getAttrOrDefault<double>(child, _Unicode(y), 0.));
0128         pos.SetZ(dd4hep::getAttrOrDefault<double>(child, _Unicode(z), 0.));
0129     }
0130     return pos;
0131 }
0132 
0133 // main
0134 static Ref_t create_detector(Detector& desc, xml::Handle_t handle, SensitiveDetector sens)
0135 {
0136     xml::DetElement detElem = handle;
0137     std::string detName = detElem.nameStr();
0138     int detID = detElem.id();
0139     DetElement det(detName, detID);
0140     sens.setType("calorimeter");
0141     // envelope
0142     Assembly assembly(detName);
0143 
0144     // module placement
0145     xml::Component plm = detElem.child(_Unicode(placements));
0146     std::map<int, int> sectorModuleNumbers;
0147     auto addModuleNumbers = [&sectorModuleNumbers] (int sector, int nmod) {
0148         auto it = sectorModuleNumbers.find(sector);
0149         if (it != sectorModuleNumbers.end()) {
0150             it->second += nmod;
0151         } else {
0152             sectorModuleNumbers[sector] = nmod;
0153         }
0154     };
0155     int sector_id = 1;
0156     for (xml::Collection_t mod(plm, _Unicode(individuals)); mod; ++mod) {
0157         auto [sector, nmod] = add_individuals(desc, assembly, mod, sens, sector_id++);
0158         addModuleNumbers(sector, nmod);
0159     }
0160     for (xml::Collection_t arr(plm, _Unicode(array)); arr; ++arr) {
0161         auto [sector, nmod] = add_array(desc, assembly, arr, sens, sector_id++);
0162         addModuleNumbers(sector, nmod);
0163     }
0164     for (xml::Collection_t disk(plm, _Unicode(disk)); disk; ++disk) {
0165         auto [sector, nmod] = add_disk(desc, assembly, disk, sens, sector_id++);
0166         addModuleNumbers(sector, nmod);
0167     }
0168     for (xml::Collection_t lines(plm, _Unicode(lines)); lines; ++lines) {
0169         auto [sector, nmod] = add_lines(desc, assembly, lines, sens, sector_id++);
0170         addModuleNumbers(sector, nmod);
0171     }
0172 
0173     for (auto [sector, nmods] : sectorModuleNumbers) {
0174         desc.add(Constant(Form((detName + "_NModules_Sector%d").c_str(), sector), std::to_string(nmods)));
0175     }
0176 
0177     // detector position and rotation
0178     auto pos = get_xml_xyz(detElem, _Unicode(position));
0179     auto rot = get_xml_xyz(detElem, _Unicode(rotation));
0180     Volume motherVol = desc.pickMotherVolume(det);
0181     Transform3D tr = Translation3D(pos.x(), pos.y(), pos.z()) * RotationZYX(rot.z(), rot.y(), rot.x());
0182     PlacedVolume envPV = motherVol.placeVolume(assembly, tr);
0183     envPV.addPhysVolID("system", detID);
0184     det.setPlacement(envPV);
0185     return det;
0186 }
0187 
0188 // helper function to build module with or w/o wrapper
0189 std::tuple<Volume, Position> build_module(Detector &desc, xml::Collection_t &plm, SensitiveDetector &sens)
0190 {
0191     auto mod = plm.child(_Unicode(module));
0192     auto sx = mod.attr<double>(_Unicode(sizex));
0193     auto sy = mod.attr<double>(_Unicode(sizey));
0194     auto sz = mod.attr<double>(_Unicode(sizez));
0195     Box modShape(sx/2., sy/2., sz/2.);
0196     auto modMat = desc.material(mod.attr<std::string>(_Unicode(material)));
0197     Volume modVol("module_vol", modShape, modMat);
0198     modVol.setSensitiveDetector(sens);
0199     modVol.setVisAttributes(desc.visAttributes(mod.attr<std::string>(_Unicode(vis))));
0200 
0201     // no wrapper
0202     if (!plm.hasChild(_Unicode(wrapper))) {
0203         return std::make_tuple(modVol, Position{sx, sy, sz});
0204     // build wrapper
0205     } else {
0206         auto wrp = plm.child(_Unicode(wrapper));
0207         auto thickness = wrp.attr<double>(_Unicode(thickness));
0208         if (thickness < 1e-12*mm) {
0209             return std::make_tuple(modVol, Position{sx, sy, sz});
0210         }
0211         auto wrpMat = desc.material(wrp.attr<std::string>(_Unicode(material)));
0212         Box wrpShape((sx + thickness)/2., (sy + thickness)/2., sz/2.);
0213         Volume wrpVol("wrapper_vol", wrpShape, wrpMat);
0214         wrpVol.placeVolume(modVol, Position(0., 0., 0.));
0215         wrpVol.setVisAttributes(desc.visAttributes(wrp.attr<std::string>(_Unicode(vis))));
0216         return std::make_tuple(wrpVol, Position{sx + thickness, sy + thickness, sz});
0217     }
0218 }
0219 
0220 // place modules, id must be provided
0221 static std::tuple<int, int> add_individuals(Detector& desc, Assembly &env, xml::Collection_t &plm, SensitiveDetector &sens, int sid)
0222 {
0223     auto [modVol, modSize] = build_module(desc, plm, sens);
0224     int sector_id = dd4hep::getAttrOrDefault<int>(plm, _Unicode(sector), sid);
0225     int nmodules = 0;
0226     for (xml::Collection_t pl(plm, _Unicode(placement)); pl; ++pl) {
0227         Position pos(dd4hep::getAttrOrDefault<double>(pl, _Unicode(x), 0.),
0228                      dd4hep::getAttrOrDefault<double>(pl, _Unicode(y), 0.),
0229                      dd4hep::getAttrOrDefault<double>(pl, _Unicode(z), 0.));
0230         Position rot(dd4hep::getAttrOrDefault<double>(pl, _Unicode(rotx), 0.),
0231                      dd4hep::getAttrOrDefault<double>(pl, _Unicode(roty), 0.),
0232                      dd4hep::getAttrOrDefault<double>(pl, _Unicode(rotz), 0.));
0233         auto mid = pl.attr<int>(_Unicode(id));
0234         Transform3D tr = Translation3D(pos.x(), pos.y(), pos.z())
0235                        * RotationZYX(rot.z(), rot.y(), rot.x());
0236         auto modPV = env.placeVolume(modVol, tr);
0237         modPV.addPhysVolID("sector", sector_id).addPhysVolID("module", mid);
0238         nmodules ++;
0239     }
0240 
0241     return {sector_id, nmodules};
0242 }
0243 
0244 // place array of modules
0245 static std::tuple<int, int> add_array(Detector& desc, Assembly &env, xml::Collection_t &plm, SensitiveDetector &sens, int sid)
0246 {
0247     auto [modVol, modSize] = build_module(desc, plm, sens);
0248     int sector_id = dd4hep::getAttrOrDefault<int>(plm, _Unicode(sector), sid);
0249     int id_begin = dd4hep::getAttrOrDefault<int>(plm, _Unicode(id_begin), 1);
0250     int nrow = plm.attr<int>(_Unicode(nrow));
0251     int ncol = plm.attr<int>(_Unicode(ncol));
0252 
0253     // compute array position
0254     double begx = -modSize.x()*ncol/2. + modSize.x()/2.;
0255     double begy = modSize.y()*nrow/2. - modSize.y()/2.;
0256 
0257     std::vector<std::pair<int, int>> removals;
0258     // get the removal list
0259     for (xml::Collection_t rm(plm, _Unicode(removal)); rm; ++rm) {
0260         removals.push_back({rm.attr<int>(_Unicode(row)), rm.attr<int>(_Unicode(col))});
0261     }
0262 
0263     // placement to mother
0264     auto pos = get_xml_xyz(plm, _Unicode(position));
0265     auto rot = get_xml_xyz(plm, _Unicode(rotation));
0266     int nmodules = 0;
0267     for (int i = 0; i < nrow; ++i) {
0268         for (int j = 0; j < ncol; ++j) {
0269             if (std::find(removals.begin(), removals.end(), std::pair<int, int>(i, j)) != removals.end()) {
0270                 continue;
0271             }
0272             double px = begx + modSize.x()*j;
0273             double py = begy - modSize.y()*i;
0274             Transform3D tr = RotationZYX(rot.z(), rot.y(), rot.x())
0275                            * Translation3D(pos.x() + px, pos.y() + py, pos.z());
0276             auto modPV = env.placeVolume(modVol, tr);
0277             modPV.addPhysVolID("sector", sector_id).addPhysVolID("module", i*ncol + j + id_begin);
0278             nmodules ++;
0279         }
0280     }
0281     return {sector_id, nmodules};
0282 }
0283 
0284 // place disk of modules
0285 static std::tuple<int, int> add_disk(Detector& desc, Assembly &env, xml::Collection_t &plm, SensitiveDetector &sens, int sid)
0286 {
0287     auto [modVol, modSize] = build_module(desc, plm, sens);
0288     int sector_id = dd4hep::getAttrOrDefault<int>(plm, _Unicode(sector), sid);
0289     int id_begin = dd4hep::getAttrOrDefault<int>(plm, _Unicode(id_begin), 1);
0290     double rmin = plm.attr<double>(_Unicode(rmin));
0291     double rmax = plm.attr<double>(_Unicode(rmax));
0292     double phimin = dd4hep::getAttrOrDefault<double>(plm, _Unicode(phimin), 0.);
0293     double phimax = dd4hep::getAttrOrDefault<double>(plm, _Unicode(phimax), 2.*M_PI);
0294 
0295     auto points = athena::geo::fillRectangles({0., 0.}, modSize.x(), modSize.y(), rmin, rmax, phimin, phimax);
0296     // placement to mother
0297     auto pos = get_xml_xyz(plm, _Unicode(position));
0298     auto rot = get_xml_xyz(plm, _Unicode(rotation));
0299     int mid = 0;
0300     for (auto &p : points) {
0301         Transform3D tr = RotationZYX(rot.z(), rot.y(), rot.x())
0302                        * Translation3D(pos.x() + p.x(), pos.y() + p.y(), pos.z());
0303         auto modPV = env.placeVolume(modVol, tr);
0304         modPV.addPhysVolID("sector", sector_id).addPhysVolID("module", id_begin + mid++);
0305     }
0306     return {sector_id, mid};
0307 }
0308 
0309 // place lines of modules (anchor point is the 0th module of this line)
0310 static std::tuple<int, int> add_lines(Detector& desc, Assembly &env, xml::Collection_t &plm, SensitiveDetector &sens, int sid)
0311 {
0312     auto [modVol, modSize] = build_module(desc, plm, sens);
0313     int sector_id = dd4hep::getAttrOrDefault<int>(plm, _Unicode(sector), sid);
0314     int id_begin = dd4hep::getAttrOrDefault<int>(plm, _Unicode(id_begin), 1);
0315     bool mirrorx = dd4hep::getAttrOrDefault<bool>(plm, _Unicode(mirrorx), false);
0316     bool mirrory = dd4hep::getAttrOrDefault<bool>(plm, _Unicode(mirrory), false);
0317     bool mirrorz = dd4hep::getAttrOrDefault<bool>(plm, _Unicode(mirrorz), false);
0318 
0319     // line placement
0320     int mid = 0;
0321     for (xml::Collection_t pl(plm, _Unicode(line)); pl; ++pl) {
0322         Position pos(dd4hep::getAttrOrDefault<double>(pl, _Unicode(x), 0.),
0323                      dd4hep::getAttrOrDefault<double>(pl, _Unicode(y), 0.),
0324                      dd4hep::getAttrOrDefault<double>(pl, _Unicode(z), 0.));
0325         Position rot(dd4hep::getAttrOrDefault<double>(pl, _Unicode(rotx), 0.),
0326                      dd4hep::getAttrOrDefault<double>(pl, _Unicode(roty), 0.),
0327                      dd4hep::getAttrOrDefault<double>(pl, _Unicode(rotz), 0.));
0328         // determine axis
0329         std::string axis = dd4hep::getAttrOrDefault(pl, _Unicode(axis), "x");
0330         std::transform(axis.begin(), axis.end(), axis.begin(), [](char c) { return std::tolower(c); });
0331         if ((axis != "x") && (axis != "y") && (axis != "z")) {
0332             std::cerr << "HomogeneousCalorimeter Error: cannot determine axis of line " << axis
0333                       << ", abort placement of this line." << std::endl;
0334             continue;
0335         }
0336 
0337         int begin = dd4hep::getAttrOrDefault<int>(pl, _Unicode(begin), 0);
0338         int nmods = pl.attr<int>(_Unicode(nmods));
0339 
0340         std::vector<Position> trans;
0341         for (int i = 0; i < nmods; ++i) {
0342             Position tran{ (axis == "x") ? pos.x() + (begin + i)*modSize.x() : pos.x(),
0343                            (axis == "y") ? pos.y() + (begin + i)*modSize.y() : pos.y(),
0344                            (axis == "z") ? pos.z() + (begin + i)*modSize.z() : pos.z() };
0345             trans.push_back(tran);
0346         }
0347 
0348         // mirror placement
0349         if (mirrorx) {
0350             int curr_size = trans.size();
0351             for (size_t i = 0; i < curr_size; ++i) {
0352                 trans.push_back(Position{-trans[i].x(), trans[i].y(), trans[i].z()});
0353             }
0354         }
0355         if (mirrory) {
0356             int curr_size = trans.size();
0357             for (size_t i = 0; i < curr_size; ++i) {
0358                 trans.push_back(Position{trans[i].x(), -trans[i].y(), trans[i].z()});
0359             }
0360         }
0361         if (mirrorz) {
0362             int curr_size = trans.size();
0363             for (size_t i = 0; i < curr_size; ++i) {
0364                 trans.push_back(Position{trans[i].x(), trans[i].y(), -trans[i].z()});
0365             }
0366         }
0367 
0368         // place volume
0369         for (auto &p : trans) {
0370             Transform3D tr = RotationZYX(rot.z(), rot.y(), rot.x())
0371                            * Translation3D(p.x(), p.y(), p.z());
0372             auto modPV = env.placeVolume(modVol, tr);
0373             modPV.addPhysVolID("sector", sector_id).addPhysVolID("module", id_begin + mid++);
0374         }
0375     }
0376     return {sector_id, mid};
0377 }
0378 //@}
0379 DECLARE_DETELEMENT(HomogeneousCalorimeter, create_detector)
0380