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
56
57
58
59
60
61
62
63
64
65
66
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
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
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
//==========================================================================
// 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 <tuple>
#include <math.h>
#include <fmt/core.h>
using namespace dd4hep;
// main
static Ref_t create_detector(Detector& desc, xml::Handle_t handle, SensitiveDetector sens)
{
using namespace std;
using namespace fmt;
xml::DetElement detElem = handle;
std::string detName = detElem.nameStr();
int detID = detElem.id();
DetElement det(detName, detID);
sens.setType("calorimeter");
auto glass_material = desc.material("PbGlass");
auto crystal_material = desc.material("PbWO4");
auto air_material = desc.material("Air");
double ROut = desc.constantAsDouble("EcalEndcapN_rmax");
double RIn = desc.constantAsDouble("EcalEndcapN_rmin");
double SizeZ = desc.constantAsDouble("EcalEndcapN_thickness");
double glass_shift_z = desc.constantAsDouble("GlassModule_z0");
double Thickness = desc.constantAsDouble("EcalEndcapN_thickness");
double trans_radius = desc.constantAsDouble("EcalEndcapNCrystal_rmax");
double Glass_ShiftZ = desc.constantAsDouble("GlassModule_z0");
double Glass_Width = desc.constantAsDouble("GlassModule_width");
double Glass_Thickness = desc.constantAsDouble("GlassModule_length");
double Glass_Gap = desc.constantAsDouble("GlassModule_wrap");
double glass_distance = desc.constantAsDouble("GlassModule_distance");
double Crystal_Width = desc.constantAsDouble("CrystalModule_width");
double Crystal_Thickness = desc.constantAsDouble("CrystalModule_length");
double Crystal_Gap = desc.constantAsDouble("CrystalModule_wrap");
double crystal_distance = desc.constantAsDouble("CrystalModule_distance");
double Crystal_shift_z = desc.constantAsDouble("CrystalModule_z0");
// RIn and ROut will define outer tube embedding the calorimeter
// centers_rin/out define the maximum radius of module centers
// so that modules are not overlapping with mother tube volume
double hypotenuse = sqrt(0.5 * glass_distance * glass_distance);
double centers_rin = RIn + hypotenuse + 1*mm;
double centers_rout = ROut - hypotenuse - 1*mm;
// envelope
// Assembly assembly(detName);
Tube outer_tube(RIn, ROut, SizeZ / 2.0, 0., 360.0 * deg);
Volume ecal_vol("negative_ecal", outer_tube, air_material);
ecal_vol.setVisAttributes(desc.visAttributes("HybridEcalOuterVis"));
double Glass_OuterR = ROut - 1 * cm ;
double Glass_InnerR = trans_radius;
glass_shift_z = Thickness / 2. - Glass_Thickness / 2.;
// Geometry of modules
Box glass_box("glass_box", Glass_Width * 0.5, Glass_Width * 0.5, Glass_Thickness * 0.5);
Volume glass_module("glass_module", glass_box, glass_material);
glass_module.setVisAttributes(desc.visAttributes("EcalEndcapNModuleVis"));
Box crystal_box("crystal_box", Crystal_Width* 0.5, Crystal_Width * 0.5, Crystal_Thickness * 0.5);
Volume crystal_module("crystal_module", crystal_box, crystal_material);
crystal_module.setVisAttributes(desc.visAttributes("EcalEndcapNModuleVis"));
// GLASS
double diameter = 2 * Glass_OuterR;
// How many towers do we have per row/columnt?
// Add a gap + diameter as if we have N towers, we have N-1 gaps;
int towersInRow = std::ceil((diameter + Glass_Gap) / (Glass_Width + Glass_Gap));
// Is it odd or even number of towersInRow
double leftTowerPos, topTowerPos;
if(towersInRow%2) {
// |
// [ ][ ][ ][ ][ ]
// ^ |
int towersInHalfRow = std::ceil(towersInRow/2.0);
topTowerPos = leftTowerPos = -towersInHalfRow * (Glass_Width + Glass_Gap);
} else {
// |
// [ ][ ][ ][ ][ ][ ]
// ^ |
int towersInHalfRow = towersInRow/2;
topTowerPos = leftTowerPos = -(towersInHalfRow - 0.5) * (Glass_Width + Glass_Gap);
}
int moduleIndex = 0;
// fmt::print("\nCE EMCAL GLASS SQUARE START\n");
// fmt::print("Glass_Thickness = {} cm;\n", Glass_Thickness / cm);
// fmt::print("Glass_Width = {} cm;\n", Glass_Width / cm);
// fmt::print("Glass_Gap = {} cm;\n", Glass_Gap / cm);
// fmt::print("Glass_InnerR = {} cm;\n", Glass_InnerR / cm);
// fmt::print("Glass_OuterR = {} cm;\n", Glass_OuterR / cm);
// fmt::print("Glass_PosZ = {} cm;\n", glass_shift_z / cm);
// fmt::print("Towers in Row/Col = {} cm;\n", glass_shift_z / cm);
// fmt::print("Top left tower pos = {:<10} {:<10} cm;\n", -leftTowerPos / cm, topTowerPos / cm);
// fmt::print("#Towers info:\n");
// fmt::print("#{:<5} {:<6} {:<3} {:<3} {:>10} {:>10} {}\n", "idx", "code", "col", "row", "x", "y", "name");
// We first do a "dry run", not really placing modules,
// but figuring out the ID scheme, number of modules, etc.
int glassModuleCount = 0;
int crystalModuleCount = 0;
int firstCrystRow = 1000000; // The first row, where crystals are started
int firstCrystCol = 1000000; // The fist column, where crystals are started
for(int rowIndex=0; rowIndex < towersInRow; rowIndex++) {
for(int colIndex=0; colIndex < towersInRow; colIndex++) {
double glass_x = leftTowerPos + colIndex * glass_distance;
double glass_y = topTowerPos + rowIndex * glass_distance;
double r = sqrt(glass_x * glass_x + glass_y * glass_y);
if (r < centers_rout && r > centers_rin) {
// we are placing something
if(r<trans_radius) {
// 4 crystal modules will be placed
crystalModuleCount+=4;
// Finding the first col and row where crystals start
// is the same algorithm as finding a minimum in array
if(colIndex<firstCrystCol) {
firstCrystCol = colIndex;
}
if(rowIndex<firstCrystRow) {
firstCrystRow = rowIndex;
}
}
else
{
// glass module will be places
glassModuleCount++;
}
}
}
}
// fmt::print("#Towers info:\n");
// fmt::print("#{:<5} {:<6} {:<3} {:<3} {:>10} {:>10} {}\n", "idx", "code", "col", "row", "x", "y", "name");
int glass_module_index = 0;
int cryst_module_index = 0;
for(int rowIndex=0; rowIndex < towersInRow; rowIndex++) {
for(int colIndex=0; colIndex < towersInRow; colIndex++) {
double glass_x = leftTowerPos + colIndex * (Glass_Width + Glass_Gap);
double glass_y = topTowerPos + rowIndex * (Glass_Width + Glass_Gap);
double r = sqrt(glass_x * glass_x + glass_y * glass_y);
if (r < centers_rout && r > centers_rin) {
int code = 1000 * rowIndex + colIndex;
std::string name = fmt::format("ce_EMCAL_glass_phys_{}", code);
if(r<trans_radius) {
// first crystal module
double crystal_x = glass_x - crystal_distance / 2;
double crystal_y = glass_y - crystal_distance / 2;
auto placement = ecal_vol.placeVolume(crystal_module, Position(crystal_x, crystal_y, Crystal_shift_z));
placement.addPhysVolID("sector", 1);
placement.addPhysVolID("module", cryst_module_index++);
// second crystal module
crystal_x = glass_x + crystal_distance / 2;
crystal_y = glass_y - crystal_distance / 2;
placement = ecal_vol.placeVolume(crystal_module, Position(crystal_x, crystal_y, Crystal_shift_z));
placement.addPhysVolID("sector", 1);
placement.addPhysVolID("module", cryst_module_index++);
// third crystal module
crystal_x = glass_x - crystal_distance / 2;
crystal_y = glass_y + crystal_distance / 2;
placement = ecal_vol.placeVolume(crystal_module, Position(crystal_x, crystal_y, Crystal_shift_z));
placement.addPhysVolID("sector", 1);
placement.addPhysVolID("module", cryst_module_index++);
// forth crystal module
crystal_x = glass_x + crystal_distance / 2;
crystal_y = glass_y + crystal_distance / 2;
placement = ecal_vol.placeVolume(crystal_module, Position(crystal_x, crystal_y, Crystal_shift_z));
placement.addPhysVolID("sector", 1);
placement.addPhysVolID("module", cryst_module_index++);
}
else
{
// glass module
auto placement = ecal_vol.placeVolume(glass_module, Position(glass_x, glass_y, glass_shift_z));
placement.addPhysVolID("sector", 2);
placement.addPhysVolID("module", glass_module_index++);
}
// fmt::print(" {:<5} {:<6} {:<3} {:<3} {:>10.4f} {:>10.4f} {}\n", towerIndex, code, colIndex, rowIndex, x / cm, y / cm, name);
glass_module_index++;
}
}
}
// fmt::print("Total Glass modules: {}\n", towerIndex);
// fmt::print("CE EMCAL GLASS END\n\n");
// detector position and rotation
auto pos = detElem.position();
auto rot = detElem.rotation();
Volume motherVol = desc.pickMotherVolume(det);
Transform3D tr(RotationZYX(rot.z(), rot.y(), rot.x()), Position(pos.x(), pos.y(), pos.z()));
PlacedVolume envPV = motherVol.placeVolume(ecal_vol, tr);
envPV.addPhysVolID("system", detID);
det.setPlacement(envPV);
return det;
}
//@}
DECLARE_DETELEMENT(HybridCalorimeter, create_detector)