File indexing completed on 2024-11-15 08:59:19
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
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
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047
0048
0049
0050
0051
0052
0053
0054
0055
0056
0057
0058
0059
0060
0061
0062
0063
0064
0065
0066
0067
0068
0069
0070
0071
0072
0073
0074
0075
0076
0077
0078
0079
0080
0081
0082
0083
0084
0085
0086
0087
0088
0089
0090
0091
0092
0093
0094
0095
0096
0097
0098
0099
0100
0101
0102
0103
0104
0105
0106
0107
0108
0109
0110
0111
0112
0113
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
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
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
0142 Assembly assembly(detName);
0143
0144
0145 xml::Component plm = detElem.child(_Unicode(placements));
0146 std::map<int, int> sectorModuleNumbers;
0147 auto addModuleNumbers = [§orModuleNumbers] (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
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
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
0202 if (!plm.hasChild(_Unicode(wrapper))) {
0203 return std::make_tuple(modVol, Position{sx, sy, sz});
0204
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
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
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
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
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
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
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
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
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
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
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
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
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