Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
//==========================================================================
// Scintillating fiber calorimeter with tower shape blocks
// reference: https://github.com/adamjaro/lmon/blob/master/calo/src/WScFiZXv3.cxx
// Support disk placement
//--------------------------------------------------------------------------
// Author: Chao Peng (ANL)
// Date: 07/19/2021
//==========================================================================
#include "GeometryHelpers.h"
#include "DD4hep/DetFactoryHelper.h"
#include <XML/Helper.h>
#include <iostream>
#include <algorithm>
#include <tuple>
#include <math.h>
using namespace dd4hep;
using Point = ROOT::Math::XYPoint;
std::tuple<Volume, Position> build_module(const Detector &desc, const xml::Component &mod_x, SensitiveDetector &sens);
// helper function to get x, y, z if defined in a xml component
template<class XmlComp>
Position get_xml_xyz(const 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)
{
xml::DetElement detElem = handle;
std::string detName = detElem.nameStr();
int detID = detElem.id();
DetElement det(detName, detID);
sens.setType("calorimeter");
auto dim = detElem.dimensions();
auto rmin = dim.rmin();
auto rmax = dim.rmax();
auto length = dim.length();
auto phimin = dd4hep::getAttrOrDefault<double>(dim, _Unicode(phimin), 0.);
auto phimax = dd4hep::getAttrOrDefault<double>(dim, _Unicode(phimax), 2.*M_PI);
// envelope
Tube envShape(rmin, rmax, length/2., phimin, phimax);
Volume env(detName + "_envelope", envShape, desc.material("Air"));
env.setVisAttributes(desc.visAttributes(detElem.visStr()));
auto [modVol, modSize] = build_module(desc, detElem.child(_Unicode(module)), sens);
double modSizeR = std::sqrt(modSize.x() * modSize.x() + modSize.y() * modSize.y());
double assembly_rwidth = modSizeR*2.;
int nas = int((rmax - rmin) / assembly_rwidth) + 1;
std::vector<Assembly> assemblies;
// calorimeter block z-offsets (as blocks are shorter than the volume length)
const double block_offset = -0.5*(length - modSize.z());
for (int i = 0; i < nas; ++i) {
Assembly assembly(detName + Form("_ring%d", i + 1));
auto assemblyPV = env.placeVolume(assembly, Position{0., 0., block_offset});
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
assemblyPV.addPhysVolID("ring", i + 1);
assemblies.emplace_back(std::move(assembly));
}
// std::cout << assemblies.size() << std::endl;
int modid = 1;
for (int ix = 0; ix < int(2.*rmax / modSize.x()) + 1; ++ix) {
double mx = modSize.x() * ix - rmax;
for (int iy = 0; iy < int(2.*rmax / modSize.y()) + 1; ++iy) {
double my = modSize.y() * iy - rmax;
double mr = std::sqrt(mx*mx + my*my);
if (mr - modSizeR >= rmin && mr + modSizeR <= rmax) {
int ias = int((mr - rmin) / assembly_rwidth);
auto &assembly = assemblies[ias];
auto modPV = assembly.placeVolume(modVol, Position(mx, my, 0.));
modPV.addPhysVolID("module", modid++);
}
}
}
desc.add(Constant(detName + "_NModules", std::to_string(modid - 1)));
for (auto &assembly : assemblies) {
assembly.ptr()->Voxelize("");
}
// 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(env, tr);
envPV.addPhysVolID("system", detID);
det.setPlacement(envPV);
return det;
}
// helper function to build module with scintillating fibers
std::tuple<Volume, Position> build_module(const Detector &desc, const xml::Component &mod_x, SensitiveDetector &sens)
{
auto sx = mod_x.attr<double>(_Unicode(sizex));
auto sy = mod_x.attr<double>(_Unicode(sizey));
auto sz = mod_x.attr<double>(_Unicode(sizez));
Box modShape(sx/2., sy/2., sz/2.);
auto modMat = desc.material(mod_x.attr<std::string>(_Unicode(material)));
Volume modVol("module_vol", modShape, modMat);
if (mod_x.hasAttr(_Unicode(vis))) {
modVol.setVisAttributes(desc.visAttributes(mod_x.attr<std::string>(_Unicode(vis))));
}
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
if (mod_x.hasChild("fiber")) {
auto fiber_x = mod_x.child(_Unicode(fiber));
auto fr = fiber_x.attr<double>(_Unicode(radius));
auto fsx = fiber_x.attr<double>(_Unicode(spacex));
auto fsy = fiber_x.attr<double>(_Unicode(spacey));
auto foff = dd4hep::getAttrOrDefault<double>(fiber_x, _Unicode(offset), 0.5*mm);
auto fiberMat = desc.material(fiber_x.attr<std::string>(_Unicode(material)));
Tube fiberShape(0., fr, sz/2.);
Volume fiberVol("fiber_vol", fiberShape, fiberMat);
fiberVol.setSensitiveDetector(sens);
// Fibers are placed in a honeycomb with the radius = sqrt(3)/2. * hexagon side length
// So each fiber is fully contained in a regular hexagon, which are placed as
// ______________________________________
// | ____ ____ |
// even: | / \ / \ |
// | ____/ \____/ \____ |
// | / \ / \ / \ |
// odd: | / \____/ \____/ \ |
// | \ / \ / \ / |
// | \____/ \____/ \____/ |
// even: | \ / \ / |
// | \____/ \____/ ___|___
// |____________________________________|___offset
// | |
// |offset
// the parameters space x and space y are used to add additional spaces between the hexagons
double fside = 2. / std::sqrt(3.) * fr;
double fdistx = 2. * fside + fsx;
double fdisty = 2. * fr + fsy;
// maximum numbers of the fibers, help narrow the loop range
int nx = int(sx / (2.*fr)) + 1;
int ny = int(sy / (2.*fr)) + 1;
// std::cout << sx << ", " << sy << ", " << fr << ", " << nx << ", " << ny << std::endl;
// place the fibers
double y0 = (foff + fside);
int nfibers = 0;
for (int iy = 0; iy < ny; ++iy) {
double y = y0 + fdisty * iy;
// about to touch the boundary
if ((sy - y) < y0) { break; }
double x0 = (iy % 2) ? (foff + fside) : (foff + fside + fdistx / 2.);
for (int ix = 0; ix < nx; ++ix) {
double x = x0 + fdistx * ix;
// about to touch the boundary
if ((sx - x) < x0) { break; }
auto fiberPV = modVol.placeVolume(fiberVol, nfibers++, Position{x - sx/2., y - sy/2., 0});
//std::cout << "(" << ix << ", " << iy << ", " << x - sx/2. << ", " << y - sy/2. << ", " << fr << "),\n";
fiberPV.addPhysVolID("fiber_x", ix + 1).addPhysVolID("fiber_y", iy + 1);
}
}