Newer
Older
//==========================================================================
// A general implementation for homogeneous calorimeter
// it supports three types of placements
// 1. Individual module placement with module dimensions and positions
// 2. Array placement with module dimensions and numbers of row and columns
// 3. Disk placement with module dimensions and (Rmin, Rmax), and (Phimin, Phimax)
// 4. Lines placement with module dimensions and (mirrorx, mirrory)
// (NOTE: anchor point is the 0th block of the line instead of line center)
//--------------------------------------------------------------------------
// Author: Chao Peng (ANL)
// Date: 06/09/2021
//==========================================================================
#include "GeometryHelpers.h"
#include "DD4hep/DetFactoryHelper.h"
#include <XML/Helper.h>
#include <iostream>
#include <algorithm>
#include <math.h>
using namespace dd4hep;
using namespace dd4hep::detail;
/** \addtogroup calorimeters Calorimeters
*/
/** \addtogroup Homogeneous Calorimeter
* \brief Type: **HomogeneousCalorimeter**.
* \author C. Peng
* \ingroup calorimeters
*
*
* \code
* <detector id="1" name="HyCal" type="HomogeneousCalorimeter" readout="EcalHits">
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
* <position x="0" y="0" z="0"/>
* <rotation x="0" y="0" z="0"/>
* <placements>
* <array nrow="34" ncol="34" sector="1">
* <position x="0" y="0" z="-9.73*cm"/>
* <module sizex="2.05*cm" sizey="2.05*cm" sizez="18*cm" vis="GreenVis" material="PbWO4"/>
* <wrapper thickness="0.015*cm" material="Epoxy" vis="WhiteVis"/>
* <removal row="16" col="16"/>
* <removal row="16" col="17"/>
* <removal row="17" col="16"/>
* <removal row="17" col="17"/>
* </array>
* <array nrow="6" ncol="24" sector="2">
* <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"/>
* <module sizex="3.8*cm" sizey="3.8*cm" sizez="45*cm" vis="BlueVis" material="PbGlass"/>
* <wrapper thickness="0.015*cm" material="Epoxy" vis="WhiteVis"/>
* </array>
* <array nrow="24" ncol="6" sector="3">
* <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"/>
* <module sizex="3.8*cm" sizey="3.8*cm" sizez="45*cm" vis="BlueVis" material="PbGlass"/>
* <wrapper thickness="0.015*cm" material="Epoxy" vis="WhiteVis"/>
* </array>
* <array nrow="6" ncol="24" sector="4">
* <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"/>
* <module sizex="3.8*cm" sizey="3.8*cm" sizez="45*cm" vis="BlueVis" material="PbGlass"/>
* <wrapper thickness="0.015*cm" material="Epoxy" vis="WhiteVis"/>
* </array>
* <array nrow="24" ncol="6" sector="3">
* <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"/>
* <module sizex="3.8*cm" sizey="3.8*cm" sizez="45*cm" vis="BlueVis" material="PbGlass"/>
* <wrapper thickness="0.015*cm" material="Epoxy" vis="WhiteVis"/>
* </array>
* </placements>
* </detector>
*
* <detector id="2" name="SomeBlocks" type="HomogeneousCalorimeter" readout="EcalHits">
* <position x="0" y="0" z="30*cm"/>
* <rotation x="0" y="0" z="0"/>
* <placements>
* <individuals sector="1"/>
* <module sizex="2.05*cm" sizey="2.05*cm" sizez="20*cm" vis="GreenVis" material="PbWO4"/>
* <wrapper thickness="0.015*cm" material="Epoxy" vis="WhiteVis"/>
* <placement x="1*cm" y="1*cm" z="0" id="1"/>
* <placement x="-1*cm" y="1*cm" z="0" id="2"/>
* <placement x="1*cm" y="-1*cm" z="0" id="3"/>
* <placement x="-1*cm" y="-1*cm" z="0" id="4"/>
* </individuals>
* </placements>
* </detector>
*
* <detector id="2" name="DiskShapeCalorimeter" type="HomogeneousCalorimeter" readout="EcalHits">
* <position x="0" y="0" z="-30*cm"/>
* <rotation x="0" y="0" z="0"/>
* <placements>
* <disk rmin="25*cm" rmax="125*cm" length="20.5*cm" phimin="0" phimax="360*degree" sector="1"/>
* <module sizex="2.05*cm" sizey="2.05*cm" sizez="20*cm" vis="GreenVis" material="PbWO4"/>
* <wrapper thickness="0.015*cm" material="Epoxy" vis="WhiteVis"/>
* </placements>
* </detector>
*
* <detector id="3" name="SomeLines" type="HomogeneousCalorimeter" readout="EcalHits">
* <position x="0" y="0" z="60*cm"/>
* <rotation x="0" y="0" z="0"/>
* <placements>
* <lines sector="1" mirrorx="true" mirrory="true"/>
* <module sizex="2.05*cm" sizey="2.05*cm" sizez="20*cm" vis="GreenVis" material="PbWO4"/>
* <wrapper thickness="0.015*cm" material="Epoxy" vis="WhiteVis"/>
* <line x="10.25*mm" y="10.25*mm" begin="8" nmods="16"/>
* <line x="10.25*mm" y="30.75*mm" begin="8" nmods="16"/>
* <line x="10.25*mm" y="51.25*mm" begin="8" nmods="16"/>
* </individuals>
* </placements>
* </detector>
* \endcode
*
* @{
*/
// headers
static void add_individuals(Detector& desc, Assembly &env, xml::Collection_t &plm, SensitiveDetector &sens, int id);
static void add_array(Detector& desc, Assembly &env, xml::Collection_t &plm, SensitiveDetector &sens, int id);
static void add_disk(Detector& desc, Assembly &env, xml::Collection_t &plm, SensitiveDetector &sens, int id);
static void add_lines(Detector& desc, Assembly &env, xml::Collection_t &plm, SensitiveDetector &sens, int id);
// helper function to get x, y, z if defined in a xml component
template<class XmlComp>
Position get_xml_xyz(XmlComp &comp, dd4hep::xml::Strng_t name)
{
Position pos(0., 0., 0.);
if (comp.hasChild(name)) {
auto child = comp.child(name);
pos.SetX(dd4hep::getAttrOrDefault<double>(child, _Unicode(x), 0.));
pos.SetY(dd4hep::getAttrOrDefault<double>(child, _Unicode(y), 0.));
pos.SetZ(dd4hep::getAttrOrDefault<double>(child, _Unicode(z), 0.));
}
return pos;
}
// main
static Ref_t create_detector(Detector& desc, xml::Handle_t handle, SensitiveDetector sens)
{
static const std::string func = "HomogeneousCalorimeter";
xml::DetElement detElem = handle;
std::string detName = detElem.nameStr();
int detID = detElem.id();
DetElement det(detName, detID);
sens.setType("calorimeter");
// envelope
Assembly assembly(detName);
// module placement
xml::Component plm = detElem.child(_Unicode(placements));
int sector = 1;
for (xml::Collection_t mod(plm, _Unicode(individuals)); mod; ++mod) {
add_individuals(desc, assembly, mod, sens, sector++);
for (xml::Collection_t arr(plm, _Unicode(array)); arr; ++arr) {
add_array(desc, assembly, arr, sens, sector++);
}
for (xml::Collection_t disk(plm, _Unicode(disk)); disk; ++disk) {
add_disk(desc, assembly, disk, sens, sector++);
}
for (xml::Collection_t lines(plm, _Unicode(lines)); lines; ++lines) {
add_lines(desc, assembly, lines, sens, sector++);
}
// detector position and rotation
auto pos = get_xml_xyz(detElem, _Unicode(position));
auto rot = get_xml_xyz(detElem, _Unicode(rotation));
Volume motherVol = desc.pickMotherVolume(det);
Transform3D tr = Translation3D(pos.x(), pos.y(), pos.z()) * RotationZYX(rot.z(), rot.y(), rot.x());
PlacedVolume envPV = motherVol.placeVolume(assembly, tr);
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
envPV.addPhysVolID("system", detID);
det.setPlacement(envPV);
return det;
}
// helper function to build module with or w/o wrapper
Volume build_module(Detector &desc, xml::Collection_t &plm, SensitiveDetector &sens, Position &dim)
{
auto mod = plm.child(_Unicode(module));
auto sx = mod.attr<double>(_Unicode(sizex));
auto sy = mod.attr<double>(_Unicode(sizey));
auto sz = mod.attr<double>(_Unicode(sizez));
dim = Position{sx, sy, sz};
Box modShape(sx/2., sy/2., sz/2.);
auto modMat = desc.material(mod.attr<std::string>(_Unicode(material)));
Volume modVol("module_vol", modShape, modMat);
modVol.setSensitiveDetector(sens);
modVol.setVisAttributes(desc.visAttributes(mod.attr<std::string>(_Unicode(vis))));
// no wrapper
if (!plm.hasChild(_Unicode(wrapper))) {
return modVol;
// build wrapper
} else {
auto wrp = plm.child(_Unicode(wrapper));
auto thickness = wrp.attr<double>(_Unicode(thickness));
if (thickness == 0.) {
return modVol;
}
auto wrpMat = desc.material(wrp.attr<std::string>(_Unicode(material)));
Box wrpShape((sx + thickness)/2., (sy + thickness)/2., sz/2.);
Volume wrpVol("wrapper_vol", wrpShape, wrpMat);
wrpVol.placeVolume(modVol, Position(0., 0., 0.));
wrpVol.setVisAttributes(desc.visAttributes(wrp.attr<std::string>(_Unicode(vis))));
dim = Position{sx + thickness, sy + thickness, sz};
return wrpVol;
}
}
// place modules, id must be provided
static void add_individuals(Detector& desc, Assembly &env, xml::Collection_t &plm, SensitiveDetector &sens, int sid)
{
Position modSize;
auto modVol = build_module(desc, plm, sens, modSize);
int sector_id = dd4hep::getAttrOrDefault<int>(plm, _Unicode(sector), sid);
for (xml::Collection_t pl(plm, _Unicode(placement)); pl; ++pl) {
Position pos(dd4hep::getAttrOrDefault<double>(pl, _Unicode(x), 0.),
dd4hep::getAttrOrDefault<double>(pl, _Unicode(y), 0.),
dd4hep::getAttrOrDefault<double>(pl, _Unicode(z), 0.));
Position rot(dd4hep::getAttrOrDefault<double>(pl, _Unicode(rotx), 0.),
dd4hep::getAttrOrDefault<double>(pl, _Unicode(roty), 0.),
dd4hep::getAttrOrDefault<double>(pl, _Unicode(rotz), 0.));
auto mid = pl.attr<int>(_Unicode(id));
Transform3D tr = Translation3D(pos.x(), pos.y(), pos.z())
* RotationZYX(rot.z(), rot.y(), rot.x());
auto modPV = env.placeVolume(modVol, tr);
modPV.addPhysVolID("sector", sector_id).addPhysVolID("module", mid);
}
}
static void add_array(Detector& desc, Assembly &env, xml::Collection_t &plm, SensitiveDetector &sens, int sid)
{
Position modSize;
auto modVol = build_module(desc, plm, sens, modSize);
int sector_id = dd4hep::getAttrOrDefault<int>(plm, _Unicode(sector), sid);
int id_begin = dd4hep::getAttrOrDefault<int>(plm, _Unicode(id_begin), 1);
int nrow = plm.attr<int>(_Unicode(nrow));
int ncol = plm.attr<int>(_Unicode(ncol));
// compute array position
double begx = -modSize.x()*ncol/2. + modSize.x()/2.;
double begy = modSize.y()*nrow/2. - modSize.y()/2.;
std::vector<std::pair<int, int>> removals;
// get the removal list
for (xml::Collection_t rm(plm, _Unicode(removal)); rm; ++rm) {
removals.push_back({rm.attr<int>(_Unicode(row)), rm.attr<int>(_Unicode(col))});
}
// placement to mother
auto pos = get_xml_xyz(plm, _Unicode(position));
auto rot = get_xml_xyz(plm, _Unicode(rotation));
for (int i = 0; i < nrow; ++i) {
for (int j = 0; j < ncol; ++j) {
if (std::find(removals.begin(), removals.end(), std::pair<int, int>(i, j)) != removals.end()) {
continue;
}
double px = begx + modSize.x()*j;
double py = begy - modSize.y()*i;
Transform3D tr = RotationZYX(rot.z(), rot.y(), rot.x())
* Translation3D(pos.x() + px, pos.y() + py, pos.z());
auto modPV = env.placeVolume(modVol, tr);
modPV.addPhysVolID("sector", sector_id).addPhysVolID("module", i*ncol + j + id_begin);
}
}
}
// place disk of modules
static void add_disk(Detector& desc, Assembly &env, xml::Collection_t &plm, SensitiveDetector &sens, int sid)
{
Position modSize;
auto modVol = build_module(desc, plm, sens, modSize);
int sector_id = dd4hep::getAttrOrDefault<int>(plm, _Unicode(sector), sid);
int id_begin = dd4hep::getAttrOrDefault<int>(plm, _Unicode(id_begin), 1);
double rmin = plm.attr<double>(_Unicode(rmin));
double rmax = plm.attr<double>(_Unicode(rmax));
double phimin = dd4hep::getAttrOrDefault<double>(plm, _Unicode(phimin), 0.);
double phimax = dd4hep::getAttrOrDefault<double>(plm, _Unicode(phimax), 2.*M_PI);
auto points = ref::utils::fillSquares({0., 0.}, modSize.x(), rmin, rmax, phimin, phimax);
// placement to mother
auto pos = get_xml_xyz(plm, _Unicode(position));
auto rot = get_xml_xyz(plm, _Unicode(rotation));
int mid = 0;
Transform3D tr = RotationZYX(rot.z(), rot.y(), rot.x())
* Translation3D(pos.x() + p.x(), pos.y() + p.y(), pos.z());
auto modPV = env.placeVolume(modVol, tr);
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
modPV.addPhysVolID("sector", sector_id).addPhysVolID("module", id_begin + mid++);
}
}
// place lines of modules (anchor point is the 0th module of this line)
static void add_lines(Detector& desc, Assembly &env, xml::Collection_t &plm, SensitiveDetector &sens, int sid)
{
Position modSize;
auto modVol = build_module(desc, plm, sens, modSize);
int sector_id = dd4hep::getAttrOrDefault<int>(plm, _Unicode(sector), sid);
int id_begin = dd4hep::getAttrOrDefault<int>(plm, _Unicode(id_begin), 1);
bool mirrorx = dd4hep::getAttrOrDefault<bool>(plm, _Unicode(mirrorx), false);
bool mirrory = dd4hep::getAttrOrDefault<bool>(plm, _Unicode(mirrory), false);
// line placement
int mid = 1;
for (xml::Collection_t pl(plm, _Unicode(line)); pl; ++pl) {
Position pos(dd4hep::getAttrOrDefault<double>(pl, _Unicode(x), 0.),
dd4hep::getAttrOrDefault<double>(pl, _Unicode(y), 0.),
dd4hep::getAttrOrDefault<double>(pl, _Unicode(z), 0.));
Position rot(dd4hep::getAttrOrDefault<double>(pl, _Unicode(rotx), 0.),
dd4hep::getAttrOrDefault<double>(pl, _Unicode(roty), 0.),
dd4hep::getAttrOrDefault<double>(pl, _Unicode(rotz), 0.));
int begin = dd4hep::getAttrOrDefault<int>(pl, _Unicode(begin), 0);
int nmods = pl.attr<int>(_Unicode(nmods));
std::vector<std::pair<double, double>> translations;
for (int i = 0; i < nmods; ++i) {
translations.push_back(std::pair<double, double>{pos.x() + (begin + i)*modSize.x(), pos.y()});
if (mirrorx) {
translations.push_back(std::pair<double, double>{-pos.x() - (begin + i)*modSize.x(), pos.y()});
}
if (mirrory) {
translations.push_back(std::pair<double, double>{pos.x() + (begin + i)*modSize.x(), -pos.y()});
}
if (mirrorx && mirrory) {
translations.push_back(std::pair<double, double>{-pos.x() - (begin + i)*modSize.x(), -pos.y()});
}
}
for (auto &p : translations) {
Transform3D tr = RotationZYX(rot.z(), rot.y(), rot.x())
* Translation3D(p.first, p.second, pos.z());
auto modPV = env.placeVolume(modVol, tr);
modPV.addPhysVolID("sector", sector_id).addPhysVolID("module", id_begin + mid++);
}