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
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
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
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
//==========================================================================
//Lumical Detector Construction
//--------------------------------------------------------------------------
//
// Author: Sun Xingyang , NJU
//==========================================================================
#include "DD4hep/DetFactoryHelper.h"
#include "DDRec/DetectorData.h"
#include "XML/Utilities.h"
#include "cmath"
#include "DDSegmentation/BitField64.h"
#include "DDSegmentation/TiledLayerGridXY.h"
#include "DDSegmentation/Segmentation.h"
#include "DDSegmentation/MultiSegmentation.h"
#include <vector>
#include <iostream>
#include "XML/Layering.h"
using namespace std;
using namespace dd4hep;
using namespace dd4hep::detail;
using dd4hep::Readout;
using dd4hep::Position;
using dd4hep::BUILD_ENVELOPE;
using dd4hep::Box;
using dd4hep::DetElement;
using dd4hep::Detector;
using dd4hep::IntersectionSolid;
using dd4hep::Material;
using dd4hep::PlacedVolume;
using dd4hep::Ref_t;
using dd4hep::Rotation3D;
using dd4hep::RotationZ;
using dd4hep::RotationZYX;
using dd4hep::SensitiveDetector;
using dd4hep::Transform3D;
using dd4hep::Trapezoid;
using dd4hep::Tube;
using dd4hep::Volume;
using dd4hep::_toString;
using dd4hep::rec::LayeredCalorimeterData;
static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector sens) {
std::cout << "This is the Lumical_v01:" << std::endl;
//get the define of LYSO cyrstals in xml
xml_det_t x_det = e;
string det_name = x_det.nameStr();
int loop = 0;
xml_coll_t LYSO_layer(x_det,_U(layer));
xml_coll_t second_LYSO_layer(x_det,_U(layer));
for(xml_coll_t c(x_det,_U(layer));c;c++){
if(loop == 3){LYSO_layer = c;}
else if(loop == 4){second_LYSO_layer=c;
break;}
loop++;
}
//get the define information of the first LYSO crystal
xml_coll_t k_0(LYSO_layer,_U(slice));
xml_comp_t component = k_0;
double dx = component.dx();
double dy = component.dy();
double dz = component.dz();
Material mat = description.material(component.materialStr());
int crystal_id,crystal_num = 0;
xml_dim_t pos = component.position();
double x = pos.x();
double y = pos.y();
double z = pos.z();
xml_dim_t rot = component.rotation();
double rotx = rot.x();
double roty = rot.y();
double rotz = rot.z();
DetElement cal(det_name, x_det.id());
// --- create an envelope volume and position it into the world ---------------------
Volume envelope = dd4hep::xml::createPlacedEnvelope( description, e, cal );
dd4hep::xml::setDetectorTypeFlag( e, cal ) ;
if( description.buildType() == BUILD_ENVELOPE ) return cal;
envelope.setVisAttributes(description, x_det.visStr());
xml_dim_t sd_typ = x_det.child(_U(sensitive));
sens.setType(sd_typ.typeStr());
Readout readout = sens.readout();
dd4hep::Segmentation seg = readout.segmentation();
dd4hep::DDSegmentation::BitField64 encoder = seg.decoder();
encoder.setValue(0) ;
dd4hep::DDSegmentation::MultiSegmentation* multiSeg =
dynamic_cast< dd4hep::DDSegmentation::MultiSegmentation*>( seg.segmentation() ) ;
dd4hep::DDSegmentation::TiledLayerGridXY* tileSeg = 0 ;
int sensitive_slice_number = -1 ;
try{
if( multiSeg ){
try{
// check if we have an entry for the subsegmentation to be used
xml_comp_t segxml = x_det.child( _Unicode( subsegmentation ) ) ;
std::string keyStr = segxml.attr<std::string>( _Unicode(key) ) ;
int keyVal = segxml.attr<int>( _Unicode(value) ) ;
encoder[ keyStr ] = keyVal ;
std::cout<<"keyStr:"<<keyStr<<"_keyVal:"<<keyVal<<"\n";
// if we have a multisegmentation that uses the slice as key, we need to know for the
// computation of the layer parameters in LayeredCalorimeterData::Layer below
if( keyStr == "layer"){sensitive_slice_number = keyVal;}
}
catch(const std::runtime_error &) {std::cerr << "Error: " << std::endl;}
// check if we have a TiledLayerGridXY segmentation :s
const dd4hep::DDSegmentation::TiledLayerGridXY* ts0 =
dynamic_cast<const dd4hep::DDSegmentation::TiledLayerGridXY*>( &multiSeg->subsegmentation( encoder.getValue() ) ) ;
tileSeg = const_cast<dd4hep::DDSegmentation::TiledLayerGridXY*>( ts0 ) ;
if( ! tileSeg ){ // if the current segmentation is not a tileSeg, we see if there is another one
for( auto s : multiSeg->subSegmentations() ){
const dd4hep::DDSegmentation::TiledLayerGridXY* ts =
dynamic_cast<const dd4hep::DDSegmentation::TiledLayerGridXY*>( s.segmentation ) ;
if( ts ) {
tileSeg = const_cast<dd4hep::DDSegmentation::TiledLayerGridXY*>( ts ) ;
break ;
}
}
}
} else {
tileSeg = dynamic_cast< dd4hep::DDSegmentation::TiledLayerGridXY*>( seg.segmentation() ) ;
}
}
catch (const std::exception& ex) {
std::cerr << "Exception caught: " << std::endl;
// Handle the exception gracefully, possibly by logging the error or providing fallback behavior
}
LayeredCalorimeterData* caloData = new LayeredCalorimeterData ;
for(int il=0;il<2; il++){
//used for reconstruction, so write a 1*1*2 layer cell size. No absorber or dead-meaterial.
dd4hep::rec::LayeredCalorimeterData::Layer _caloLayer;
_caloLayer.distance = 560+il*80*mm;
_caloLayer.phi0 = 0;
_caloLayer.absorberThickness = 0;
_caloLayer.inner_nRadiationLengths = 0.01;
_caloLayer.inner_nInteractionLengths = 0.01;
_caloLayer.outer_nRadiationLengths = 0.01;
_caloLayer.outer_nInteractionLengths = 0.01;
_caloLayer.inner_thickness = 3*mm; //1cm
_caloLayer.outer_thickness = 3*mm; //1cm
_caloLayer.sensitive_thickness = 2*3*mm; //2cm
_caloLayer.cellSize0 = 3*mm; //1cm
_caloLayer.cellSize1 = 3*mm; //1cm
caloData->layers.push_back(_caloLayer);
}
caloData->layoutType = LayeredCalorimeterData::BarrelLayout ;
caloData->inner_symmetry = 8 ;
caloData->outer_symmetry = 8 ;
caloData->phi0 = 0 ; // hardcoded
// extent of the calorimeter in the r-z-plane [ rmin, rmax, zmin, zmax ] in mm.
caloData->extent[0] = 10*mm ;
caloData->extent[1] = 50*mm;
caloData->extent[2] = 0. ;
caloData->extent[3] = 700*mm ;
//loop in order to build symmetry geo
for(int kz=1;kz>=-1;kz-=2){
/////////////////////////////////////////////////////////////////
//build Disk and flange
std::vector<double> cellSizeVector = seg.cellDimensions( encoder.getValue() ); //Assume uniform cell sizes, provide dummy cellID
int layer_num=1;// count the number of layers in order to define the ID of objects
for(int ky = 1;ky>=-1;ky-=2){
int logical_disk_id = 1; //count the number of disks
int l=0;//stand for the ID number of module
for(xml_coll_t c(x_det,_U(layer));c;c++){
string module_name = _toString(l+1,"_module%d")+_toString(0,"_stave%d");
int slice_number = 0;//used to count the number of slices we have placed
xml_comp_t x_layer = c;
string layer_name = det_name+module_name+_toString(layer_num,"_layer%d");
encoder["layer"] = logical_disk_id ;
cellSizeVector = seg.segmentation()->cellDimensions( encoder.getValue() );
int slice_num = 0;//used to distinguish which slice we are building
if(l>=2){
//build flange
if(ky==-1)
{
break;
}
slice_number=1;
for(xml_coll_t k(x_layer,_U(slice)); k; k++){
xml_comp_t x_slice = k;
string slice_name = layer_name+_toString(slice_num,"_slice%d");
DetElement slice(slice_name,_toString(slice_num,"_slice%d"),x_det.id());
Material slice_material = description.material(x_slice.materialStr());
Volume slice_vol(slice_name,Tube(x_slice.rmin(),x_slice.rmax(),x_slice.z()/2,x_slice.phi1(),x_slice.phi2()),slice_material);
slice_vol.setVisAttributes(description,x_slice.visStr());
PlacedVolume slice_phv = envelope.placeVolume(slice_vol,Position(0*mm,0*mm,kz*x_slice.position().z()));
slice_phv.addPhysVolID("side",kz).addPhysVolID("module",4).addPhysVolID("layer",0 ).addPhysVolID("slice",layer_num);
if(ky==1 && kz==1){
std::cout<<"Flange"<<layer_num<<":"<<"zStart = "<<x_slice.position().z()-x_slice.z()/2 << " zEnd = "<<x_slice.position().z()+x_slice.z()/2 <<" rmin ="<<x_slice.rmin()<<" rmax = "<<x_slice.rmax()<<"position of circle center = "<<x_slice.position().x()<<" "<<x_slice.position().y()<<" material = "<<x_slice.materialStr()<<"\n";
}
slice.setPlacement(slice_phv);
logical_disk_id++;
slice_number++;
}
break;
}
for(xml_coll_t k(x_layer,_U(slice)); k; k++) {
xml_comp_t x_slice = k;
string slice_name = layer_name + _toString(slice_number,"_slice%d");
Material slice_material = description.material(x_slice.materialStr());
DetElement slice(layer_name,_toString(slice_number,"slice%d"),x_det.id());
if(slice_num%2==0){
Volume slice_vol(slice_name,Box(x_slice.dx()/2,x_slice.dy()/2,x_slice.dz()/2),slice_material);
slice_vol.setSensitiveDetector(sens);
slice_vol.setVisAttributes(description,x_slice.visStr());
PlacedVolume slice_phv = envelope.placeVolume(slice_vol,Position(x_slice.position().x(),ky*x_slice.position().y(),kz*x_slice.position().z()));
slice_phv.addPhysVolID("side",kz).addPhysVolID("stave",ky).addPhysVolID("module",l+1).addPhysVolID("layer",layer_num ).addPhysVolID("slice",slice_number);
if(ky==1 && kz==1){
std::cout<<"Disk_Box:"<<"zStart = "<<x_slice.position().z()-x_slice.dz()/2 << " zEnd = "<<x_slice.position().z()+x_slice.dz()/2 <<" dy = "<<x_slice.dy()<<" dx = "<<x_slice.dx()<<" position of mass center = "<<x_slice.position().x()<<" "<<x_slice.position().y()<<" material = "<<x_slice.materialStr()<<"\n";
}
slice.setPlacement(slice_phv);
slice_num++;
slice_number++;
}
else{
for (int kx=1;kx>=-1;kx-=2){
double phi1 = 0*deg;
if(kx==1){phi1 =0*deg;}
else{phi1 = 90*deg;}
if(ky==-1){phi1+=90*deg;
if(kx==1)phi1+=180*deg;}//invert the geometry
Volume slice_vol(slice_name,Tube(x_slice.rmin(),x_slice.rmax(),x_slice.z()/2,phi1,phi1+90*deg),slice_material);
slice_vol.setSensitiveDetector(sens);
slice_vol.setVisAttributes(description,x_slice.visStr());
PlacedVolume slice_phv = envelope.placeVolume(slice_vol,Position(kx*x_slice.position().x(),ky*x_slice.position().y(),kz*x_slice.position().z()));
slice_phv.addPhysVolID("side",kz).addPhysVolID("stave",ky).addPhysVolID("module",l+1).addPhysVolID("layer",layer_num ).addPhysVolID("slice",slice_number);
slice.setPlacement(slice_phv);
if(ky==1 && kz==1){
std::cout<<"Disk_tube:"<<"zStart = "<<x_slice.position().z()-x_slice.z()/2 << " zEnd = "<<x_slice.position().z()+x_slice.z()/2 <<" r = "<<x_slice.rmax()<<"position of circle center = "<<x_slice.position().x()<<" "<<x_slice.position().y()<<" material = "<<x_slice.materialStr()<<"\n";
}
slice_number++;
}
slice_num++;
}
}
l=l+1;
logical_disk_id++;
}
layer_num++;
}
//////////////////////////////////////////////////////////////////
//build 1st LYSO
cellSizeVector = seg.cellDimensions( encoder.getValue() );
for(int ky=-1;ky<=1;ky=ky+2) {
string module_name = _toString(3,"_module%d")+_toString(ky,"_stave%d");
for (int i = 0; i < 14; i++) {
encoder["layer"] = 4 ;
cellSizeVector = seg.segmentation()->cellDimensions( encoder.getValue() );
int num = 0;//calculate the number of crystals allowed to be placed in each line
crystal_id=0;//stand for the id of LYSO crystals
double yc = 12.0 + dy * (i + 1);
double xc = dx / mm;
num = (int) 2 * ((std::sqrt(56 * 56 - yc * yc) / xc) - 1);
string layer_name = det_name+ module_name+_toString(kz*(i+1),"_layer%d");
int j = 0;
int half = 0;
if (num % 2 == 1) {
half = (num - 1) / 2; // half amount of total number each line
j = -half;
while (j <= half) {
string crystal_name = layer_name +_toString(crystal_id,"_slice%d");
DetElement crystalDE(layer_name, _toString(crystal_num,"_slice%d"), x_det.id());
Box crystalBox(0.5*dx, 0.5*dy, 0.5*dz);
Volume crystalVol(crystal_name, crystalBox, mat);
crystalVol.setVisAttributes(description, component.visStr());
crystalVol.setSensitiveDetector(sens);
Transform3D transform(RotationZYX(rotz, roty, rotx), Translation3D(x + j * dx, ky*(y - 0.5*dy-i * dy), kz*z));
PlacedVolume LYSO_phv =envelope.placeVolume(crystalVol, transform);
if (x_det.hasAttr(_U(id))) LYSO_phv.addPhysVolID("system", x_det.id());
LYSO_phv.addPhysVolID("side",kz).addPhysVolID("stave",ky).addPhysVolID("module",3).addPhysVolID("layer",i+1 ).addPhysVolID("slice",crystal_id);
if(ky==1 && kz==1){
std::cout << "LYSO_pixel:" << "zStart = " << z - dz / 2 << " zEnd = "<< z + dz / 2 <<" dy = "<< dy << " dx = "<< dx << " position of mass center = "<< x + j * dx <<" "<< ky * (y - 0.5*dy-i * dy)<<"\n";
}
crystalDE.setPlacement(LYSO_phv);
j++;
crystal_id++;
crystal_num++;
}
} else {
half = num / 2;
j = -half;
while (j < half) {
string crystal_name = layer_name+_toString(crystal_id,"_slice%d");
DetElement crystalDE(layer_name, _toString(crystal_num,"_slice%d"), x_det.id());
Box crystalBox(0.5*dx, 0.5*dy, 0.5*dz);
Volume crystalVol(crystal_name, crystalBox, mat);
crystalVol.setVisAttributes(description, component.visStr());
crystalVol.setSensitiveDetector(sens);
Transform3D transform(RotationZYX(rotz, roty, rotx),
Translation3D(x + (j + 0.5) * dx, ky*(y -0.5*dy- i * dy), kz*z));
PlacedVolume LYSO_phv = envelope.placeVolume(crystalVol, transform);
if (x_det.hasAttr(_U(id))) LYSO_phv.addPhysVolID("system", x_det.id());
LYSO_phv.addPhysVolID("side",kz).addPhysVolID("stave",ky).addPhysVolID("module",3).addPhysVolID("layer",i+1 ).addPhysVolID("slice",crystal_id);
if(ky==1 && kz==1){
std::cout<<"LYSO_pixel:"<<"zStart = "<<z-dz/2<< " zEnd = "<<z+dz/2 <<" dy = "<< dy <<" dx = "<<dx<<" position of mass center = "<<x + j * dx<<" "<<ky*(y - 0.5*dy-i * dy)<<"\n";
}
crystalDE.setPlacement(LYSO_phv);
j++;
crystal_id++;
crystal_num++;
}
}
}
}
//////////////////////////////////////////////////////////////////////
// build 2nd LYSO
xml_coll_t k_second(second_LYSO_layer,_U(slice));
xml_comp_t component_second = k_second;
double dx2 = component_second.dx();
double dy2 = component_second.dy();
double dz2 = component_second.dz();
Material mat2 = description.material(component_second.materialStr());
xml_dim_t pos_second = component_second.position();
double x2 = pos_second.x();
double y2 = pos_second.y();
double z2 = pos_second.z();
for(int ky=-1;ky<=1;ky=ky+2) {
string module_name = _toString(5,"_module%d")+_toString(ky,"_stave%d");
for (int i = 0; i < 10; i++) {
encoder["layer"] = 4 ;
cellSizeVector = seg.segmentation()->cellDimensions( encoder.getValue() );
int num = 0;//calculate the number of crystals allowed to be placed in each line
crystal_id=0;//stand for the id of LYSO crystals
double yc =12.0+ dy2/mm * (i+1);
double xc = dx2 / mm;
num = (int) 2 * ((std::sqrt(100.0*100.0 - yc * yc) / xc) - 1);
string layer_name = det_name+ module_name+_toString(kz*(i+1),"_layer%d");
int j = 0;
int half = 0;
if (num % 2 == 1) {
half = (num - 1) / 2;
j = -half;
while (j <= half) {
string crystal_name = layer_name +_toString(crystal_id,"_slice%d");
DetElement crystalDE(layer_name, _toString(crystal_num,"_slice%d"), x_det.id());
Box crystalBox(0.5*dx2, 0.5*dy2, 0.5*dz2);
Volume crystalVol(crystal_name, crystalBox, mat2);
crystalVol.setVisAttributes(description, component_second.visStr());
crystalVol.setSensitiveDetector(sens);
Transform3D transform(RotationZYX(0, 0, 0), Translation3D(x2 + j * dx2, ky*(y2 -0.5*dy2- i * dy2), kz*z2));
PlacedVolume LYSO_phv =envelope.placeVolume(crystalVol, transform);
if (x_det.hasAttr(_U(id))) LYSO_phv.addPhysVolID("system", x_det.id());
LYSO_phv.addPhysVolID("side",kz).addPhysVolID("stave",ky).addPhysVolID("module",5).addPhysVolID("layer",i+1 ).addPhysVolID("slice",crystal_id);
if(ky==1 && kz==1){
std::cout<<"LYSO_pixel:"<<"zStart = "<<z2-dz2/2<< " zEnd = "<<z2+dz2/2 <<" dy = "<< dy2 <<" dx = "<<dx2<<" position of mass center = "<<x2 + j * dx2<<" "<<ky*(y2 - 0.5*dy2-i * dy2)<<"\n";
}
crystalDE.setPlacement(LYSO_phv);
j++;
crystal_id++;
crystal_num++;
}
} else {
half = num / 2;
j = -half;
while (j < half) {
string crystal_name = layer_name+_toString(crystal_id,"_slice%d");
DetElement crystalDE(layer_name, _toString(crystal_num,"_slice%d"), x_det.id());
Box crystalBox(0.5*dx2, 0.5*dy2, 0.5*dz2);
Volume crystalVol(crystal_name, crystalBox, mat2);
crystalVol.setVisAttributes(description, component_second.visStr());
crystalVol.setSensitiveDetector(sens);
Transform3D transform(RotationZYX(0, 0, 0),
Translation3D(x2 + (j + 0.5) * dx2, ky*(y2 -0.5*dy2- i * dy2), kz*z2));
PlacedVolume LYSO_phv = envelope.placeVolume(crystalVol, transform);
if (x_det.hasAttr(_U(id))) LYSO_phv.addPhysVolID("system", x_det.id());
LYSO_phv.addPhysVolID("side",kz).addPhysVolID("stave",ky).addPhysVolID("module",5).addPhysVolID("layer",i+1 ).addPhysVolID("slice",crystal_id);
if(ky==1 && kz==1){
std::cout<<"LYSO_pixel:"<<"zStart = "<<z2-dz2/2<< " zEnd = "<<z2+dz2/2 <<" dy = "<< dy2 <<" dx = "<<dx2<<" position of mass center = "<<x2 + j * dx2<<" "<<ky*(y2 - 0.5*dy2-i * dy2)<<"\n";
}
crystalDE.setPlacement(LYSO_phv);
j++;
crystal_id++;
crystal_num++;
}
}
}
}
}
cal.addExtension< LayeredCalorimeterData >( caloData ) ;
return cal;
}
DECLARE_DETELEMENT(Lumical_v01, create_detector)