Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • EIC/detectors/athena
  • zwzhao/athena
  • FernandoTA/athena
  • palspeic/athena
4 results
Show changes
Hide whitespace changes
Inline Side-by-side
Showing
with 4171 additions and 205 deletions
src/BarrelCalorimeter_geo.cpp
View file @ 1008a184
...@@ -12,7 +12,7 @@ ...@@ -12,7 +12,7 @@
//========================================================================== //==========================================================================
// //
// Specialized generic detector constructor // Specialized generic detector constructor
// //
//========================================================================== //==========================================================================
#include "DD4hep/DetFactoryHelper.h" #include "DD4hep/DetFactoryHelper.h"
#include "XML/Layering.h" #include "XML/Layering.h"
...@@ -34,7 +34,11 @@ static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector s ...@@ -34,7 +34,11 @@ static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector s
double inner_r = x_dim.rmin(); double inner_r = x_dim.rmin();
double dphi = (2*M_PI/nsides); double dphi = (2*M_PI/nsides);
double hphi = dphi/2; double hphi = dphi/2;
double mod_z = layering.totalThickness(); double support_thickness = 0.0;
if(x_staves.hasChild("support")){
support_thickness = getAttrOrDefault(x_staves.child(_U(support)), _U(thickness), 5.0 * cm);
}
double mod_z = layering.totalThickness() + support_thickness;
double outer_r = inner_r + mod_z; double outer_r = inner_r + mod_z;
double totThick = mod_z; double totThick = mod_z;
double offset = x_det.attr<double>(_Unicode(offset)); double offset = x_det.attr<double>(_Unicode(offset));
...@@ -50,14 +54,14 @@ static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector s ...@@ -50,14 +54,14 @@ static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector s
DetElement stave_det("stave0",det_id); DetElement stave_det("stave0",det_id);
double dx = 0.0; //mod_z / std::sin(dphi); // dx per layer double dx = 0.0; //mod_z / std::sin(dphi); // dx per layer
// Compute the top and bottom face measurements. // Compute the top and bottom face measurements.
double trd_x2 = (2 * std::tan(hphi) * outer_r - dx)/2 - tolerance; double trd_x2 = (2 * std::tan(hphi) * outer_r - dx)/2 - tolerance;
double trd_x1 = (2 * std::tan(hphi) * inner_r + dx)/2 - tolerance; double trd_x1 = (2 * std::tan(hphi) * inner_r + dx)/2 - tolerance;
double trd_y1 = x_dim.z()/2 - tolerance; double trd_y1 = x_dim.z()/2 - tolerance;
double trd_y2 = trd_y1; double trd_y2 = trd_y1;
double trd_z = mod_z/2 - tolerance; double trd_z = mod_z/2 - tolerance;
// Create the trapezoid for the stave. // Create the trapezoid for the stave.
Trapezoid trd(trd_x1, // Outer side, i.e. the "short" X side. Trapezoid trd(trd_x1, // Outer side, i.e. the "short" X side.
trd_x2, // Inner side, i.e. the "long" X side. trd_x2, // Inner side, i.e. the "long" X side.
...@@ -66,6 +70,65 @@ static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector s ...@@ -66,6 +70,65 @@ static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector s
trd_z); // Thickness, in Y for top stave, when rotated. trd_z); // Thickness, in Y for top stave, when rotated.
Volume mod_vol("stave",trd,air); Volume mod_vol("stave",trd,air);
double l_pos_z = -(layering.totalThickness() / 2) - support_thickness/2.0;
//double trd_x2_support = trd_x1;
double trd_x1_support = (2 * std::tan(hphi) * outer_r - dx- support_thickness)/2 - tolerance;
Solid support_frame_s;
// optional stave support
if(x_staves.hasChild("support")){
xml_comp_t x_support = x_staves.child(_U(support));
// is the support on the inside surface?
bool is_inside_support = getAttrOrDefault<bool>(x_support, _Unicode(inside), true);
// number of "beams" running the length of the stave.
int n_beams = getAttrOrDefault<int>(x_support, _Unicode(n_beams), 3);
double beam_thickness = support_thickness / 4.0; // maybe a parameter later...
trd_x1_support = (2 * std::tan(hphi) * (outer_r - support_thickness + beam_thickness)) / 2 - tolerance;
double grid_size = getAttrOrDefault(x_support, _Unicode(grid_size), 25.0 * cm);
double beam_width = 2.0 * trd_x1_support / (n_beams + 1); // quick hack to make some gap between T beams
double cross_beam_thickness = support_thickness/4.0;
//double trd_x1_support = (2 * std::tan(hphi) * (inner_r + beam_thickness)) / 2 - tolerance;
double trd_x2_support = trd_x2;
int n_cross_supports = std::floor((trd_y1-cross_beam_thickness)/grid_size);
Box beam_vert_s(beam_thickness / 2.0 - tolerance, trd_y1, support_thickness / 2.0 - tolerance);
Box beam_hori_s(beam_width / 2.0 - tolerance, trd_y1, beam_thickness / 2.0 - tolerance);
UnionSolid T_beam_s(beam_vert_s, beam_hori_s, Position(0, 0, -support_thickness / 2.0 + beam_thickness / 2.0));
// cross supports
Trapezoid trd_support(trd_x1_support,trd_x2_support,
beam_thickness / 2.0 - tolerance, beam_thickness / 2.0 - tolerance,
support_thickness / 2.0 - tolerance - cross_beam_thickness/2.0);
UnionSolid support_array_start_s(T_beam_s,trd_support,Position(0,0,cross_beam_thickness/2.0));
for (int isup = 0; isup < n_cross_supports; isup++) {
support_array_start_s = UnionSolid(support_array_start_s, trd_support, Position(0, -1.0 * isup * grid_size, cross_beam_thickness/2.0));
support_array_start_s = UnionSolid(support_array_start_s, trd_support, Position(0, 1.0 * isup * grid_size, cross_beam_thickness/2.0));
}
support_array_start_s =
UnionSolid(support_array_start_s, beam_hori_s,
Position(-1.8 * 0.5*(trd_x1+trd_x2_support) / n_beams, 0, -support_thickness / 2.0 + beam_thickness / 2.0));
support_array_start_s =
UnionSolid(support_array_start_s, beam_hori_s,
Position(1.8 * 0.5*(trd_x1+trd_x2_support) / n_beams, 0, -support_thickness / 2.0 + beam_thickness / 2.0));
support_array_start_s =
UnionSolid(support_array_start_s, beam_vert_s, Position(-1.8 * 0.5*(trd_x1+trd_x2_support) / n_beams, 0, 0));
support_array_start_s =
UnionSolid(support_array_start_s, beam_vert_s, Position(1.8 * 0.5*(trd_x1+trd_x2_support) / n_beams, 0, 0));
support_frame_s = support_array_start_s;
Material support_mat = description.material(x_support.materialStr());
Volume support_vol("support_frame_v", support_frame_s, support_mat);
support_vol.setVisAttributes(description,x_support.visStr());
// figure out how to best place
//auto pv = mod_vol.placeVolume(support_vol, Position(0.0, 0.0, l_pos_z + support_thickness / 2.0));
auto pv = mod_vol.placeVolume(support_vol, Position(0.0, 0.0, -l_pos_z - support_thickness / 2.0));
}
//l_pos_z += support_thickness;
sens.setType("calorimeter"); sens.setType("calorimeter");
{ // ===== buildBarrelStave(description, sens, module_volume) ===== { // ===== buildBarrelStave(description, sens, module_volume) =====
...@@ -73,7 +136,6 @@ static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector s ...@@ -73,7 +136,6 @@ static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector s
double stave_z = trd_y1; double stave_z = trd_y1;
double tan_hphi = std::tan(hphi); double tan_hphi = std::tan(hphi);
double l_dim_x = trd_x1; // Starting X dimension for the layer. double l_dim_x = trd_x1; // Starting X dimension for the layer.
double l_pos_z = -(layering.totalThickness() / 2);
// Loop over the sets of layer elements in the detector. // Loop over the sets of layer elements in the detector.
int l_num = 1; int l_num = 1;
...@@ -112,10 +174,10 @@ static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector s ...@@ -112,10 +174,10 @@ static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector s
slice.setPlacement(slice_phv); slice.setPlacement(slice_phv);
// Increment Z position of slice. // Increment Z position of slice.
s_pos_z += s_thick; s_pos_z += s_thick;
// Increment slice number. // Increment slice number.
++s_num; ++s_num;
} }
// Set region, limitset, and vis of layer. // Set region, limitset, and vis of layer.
layer.setAttributes(description,l_vol,x_layer.regionStr(),x_layer.limitsStr(),x_layer.visStr()); layer.setAttributes(description,l_vol,x_layer.regionStr(),x_layer.limitsStr(),x_layer.visStr());
...@@ -126,7 +188,7 @@ static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector s ...@@ -126,7 +188,7 @@ static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector s
// Increment to next layer Z position. // Increment to next layer Z position.
double xcut = l_thickness * tan_hphi; double xcut = l_thickness * tan_hphi;
l_dim_x += xcut; l_dim_x += xcut;
l_pos_z += l_thickness; l_pos_z += l_thickness;
++l_num; ++l_num;
} }
} }
...@@ -135,6 +197,7 @@ static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector s ...@@ -135,6 +197,7 @@ static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector s
// Set stave visualization. // Set stave visualization.
if ( x_staves ) { if ( x_staves ) {
mod_vol.setVisAttributes(description.visAttributes(x_staves.visStr())); mod_vol.setVisAttributes(description.visAttributes(x_staves.visStr()));
} }
// Phi start for a stave. // Phi start for a stave.
double phi = M_PI / nsides; double phi = M_PI / nsides;
...@@ -161,5 +224,5 @@ static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector s ...@@ -161,5 +224,5 @@ static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector s
return sdet; return sdet;
} }
DECLARE_DETELEMENT(refdet_EcalBarrel,create_detector) DECLARE_DETELEMENT(athena_EcalBarrel,create_detector)
DECLARE_DETELEMENT(refdet_HcalBarrel,create_detector) DECLARE_DETELEMENT(athena_HcalBarrel,create_detector)
src/BarrelTrackerWithFrame_geo.cpp 0 → 100644
View file @ 1008a184
/** \addtogroup Trackers Trackers
* \brief Type: **BarrelTrackerWithFrame**.
* \author W. Armstrong
*
* \ingroup trackers
*
* @{
*/
#include "DD4hep/DetFactoryHelper.h"
#include "DD4hep/Printout.h"
#include "DD4hep/Shapes.h"
#include "DDRec/Surface.h"
#include "DDRec/DetectorData.h"
#include "XML/Layering.h"
#include "XML/Utilities.h"
#include <array>
#if defined(USE_ACTSDD4HEP)
#include "ActsDD4hep/ActsExtension.hpp"
#include "ActsDD4hep/ConvertMaterial.hpp"
#else
#include "Acts/Plugins/DD4hep/ActsExtension.hpp"
#include "Acts/Plugins/DD4hep/ConvertDD4hepMaterial.hpp"
#endif
using namespace std;
using namespace dd4hep;
using namespace dd4hep::rec;
using namespace dd4hep::detail;
/** Barrel Tracker with space frame.
*
* - Optional "support" tag within the detector element.
*
* The shapes are created using createShape which can be one of many basic geomtries.
* See the examples Check_shape_*.xml in
* [dd4hep's examples/ClientTests/compact](https://github.com/AIDASoft/DD4hep/tree/master/examples/ClientTests/compact)
* directory.
*
*
* - Optional "frame" tag within the module element.
*
* \ingroup trackers
*
* \code
* \endcode
*
*
* @author Whitney Armstrong
*/
static Ref_t create_BarrelTrackerWithFrame(Detector& description, xml_h e, SensitiveDetector sens) {
typedef vector<PlacedVolume> Placements;
xml_det_t x_det = e;
Material air = description.air();
int det_id = x_det.id();
string det_name = x_det.nameStr();
DetElement sdet(det_name, det_id);
map<string, Volume> volumes;
map<string, Placements> sensitives;
map<string, std::vector<VolPlane>> volplane_surfaces;
map<string, std::array<double, 2>> module_thicknesses;
PlacedVolume pv;
dd4hep::xml::Dimension dimensions(x_det.dimensions());
// ACTS extension
{
Acts::ActsExtension* detWorldExt = new Acts::ActsExtension();
detWorldExt->addType("barrel", "detector");
// Add the volume boundary material if configured
for (xml_coll_t bmat(x_det, _Unicode(boundary_material)); bmat; ++bmat) {
xml_comp_t x_boundary_material = bmat;
Acts::xmlToProtoSurfaceMaterial(x_boundary_material, *detWorldExt, "boundary_material");
}
sdet.addExtension<Acts::ActsExtension>(detWorldExt);
}
Tube topVolumeShape(dimensions.rmin(), dimensions.rmax(), dimensions.length() * 0.5);
Volume assembly(det_name,topVolumeShape,air);
sens.setType("tracker");
// Loop over the suports
for (xml_coll_t su(x_det, _U(support)); su; ++su) {
xml_comp_t x_support = su;
double support_thickness = getAttrOrDefault(x_support, _U(thickness), 2.0 * mm);
double support_length = getAttrOrDefault(x_support, _U(length), 2.0 * mm);
double support_rmin = getAttrOrDefault(x_support, _U(rmin), 2.0 * mm);
double support_zstart = getAttrOrDefault(x_support, _U(zstart), 2.0 * mm);
std::string support_name = getAttrOrDefault<std::string>(x_support, _Unicode(name), "support_tube");
std::string support_vis = getAttrOrDefault<std::string>(x_support, _Unicode(vis), "AnlRed");
xml_dim_t pos (x_support.child(_U(position), false));
xml_dim_t rot (x_support.child(_U(rotation), false));
Solid support_solid;
if(x_support.hasChild("shape")){
xml_comp_t shape(x_support.child(_U(shape)));
string shape_type = shape.typeStr();
support_solid = xml::createShape(description, shape_type, shape);
} else {
support_solid = Tube(support_rmin, support_rmin + support_thickness, support_length / 2);
}
Transform3D tr = Transform3D(Rotation3D(),Position(0,0,(support_zstart + support_length / 2)));
if ( pos.ptr() && rot.ptr() ) {
Rotation3D rot3D(RotationZYX(rot.z(0),rot.y(0),rot.x(0)));
Position pos3D(pos.x(0),pos.y(0),pos.z(0));
tr = Transform3D(rot3D, pos3D);
}
else if ( pos.ptr() ) {
tr = Transform3D(Rotation3D(),Position(pos.x(0),pos.y(0),pos.z(0)));
}
else if ( rot.ptr() ) {
Rotation3D rot3D(RotationZYX(rot.z(0),rot.y(0),rot.x(0)));
tr = Transform3D(rot3D,Position());
}
Material support_mat = description.material(x_support.materialStr());
Volume support_vol(support_name, support_solid, support_mat);
support_vol.setVisAttributes(description.visAttributes(support_vis));
pv = assembly.placeVolume(support_vol, tr);
// pv = assembly.placeVolume(support_vol, Position(0, 0, support_zstart + support_length / 2));
}
// loop over the modules
for (xml_coll_t mi(x_det, _U(module)); mi; ++mi) {
xml_comp_t x_mod = mi;
string m_nam = x_mod.nameStr();
if (volumes.find(m_nam) != volumes.end()) {
printout(ERROR, "BarrelTrackerWithFrame", string((string("Module with named ") + m_nam + string(" already exists."))).c_str() );
throw runtime_error("Logics error in building modules.");
}
int ncomponents = 0;
int sensor_number = 1;
double total_thickness = 0;
// Compute module total thickness from components
xml_coll_t ci(x_mod, _U(module_component));
for (ci.reset(), total_thickness = 0.0; ci; ++ci) {
total_thickness += xml_comp_t(ci).thickness();
}
// the module assembly volume
Assembly m_vol( m_nam );
volumes[m_nam] = m_vol;
m_vol.setVisAttributes(description.visAttributes(x_mod.visStr()));
// Optional module frame.
if(x_mod.hasChild("frame")){
xml_comp_t m_frame = x_mod.child(_U(frame));
//xmleles[m_nam] = x_mod;
double frame_thickness = m_frame.thickness();
double frame_width = m_frame.width();
double frame_height = getAttrOrDefault<double>(m_frame, _U(height), 5.0 * mm);
double tanth = frame_height/(frame_width/2.0);
double frame_height2 = frame_height-frame_thickness-frame_thickness/tanth;
double frame_width2 = 2.0*frame_height2/tanth;
Trd1 moduleframe_part1(frame_width / 2, 0.001 * mm, m_frame.length() / 2,
frame_height / 2);
Trd1 moduleframe_part2(frame_width2/2, 0.001 * mm,
m_frame.length() / 2 + 0.01 * mm, frame_height2/2);
SubtractionSolid moduleframe(moduleframe_part1, moduleframe_part2,Position(0.0,frame_thickness,0.0));
Volume v_moduleframe(m_nam+"_vol", moduleframe, description.material(m_frame.materialStr()));
v_moduleframe.setVisAttributes(description, m_frame.visStr());
m_vol.placeVolume(v_moduleframe, Position(0.0, 0.0, frame_height / 2 + total_thickness / 2.0));
}
double thickness_so_far = 0.0;
double thickness_sum = -total_thickness/2.0;
for (xml_coll_t ci(x_mod, _U(module_component)); ci; ++ci, ++ncomponents) {
xml_comp_t x_comp = ci;
xml_comp_t x_pos = x_comp.position(false);
xml_comp_t x_rot = x_comp.rotation(false);
const string c_nam = _toString(ncomponents, "component%d");
Box c_box(x_comp.width() / 2, x_comp.length() / 2, x_comp.thickness() / 2);
Volume c_vol(c_nam, c_box, description.material(x_comp.materialStr()));
// Utility variable for the relative z-offset based off the previous components
const double zoff = thickness_sum+x_comp.thickness() / 2.0;
if (x_pos && x_rot) {
Position c_pos(x_pos.x(0), x_pos.y(0), x_pos.z(0) + zoff);
RotationZYX c_rot(x_rot.z(0), x_rot.y(0), x_rot.x(0));
pv = m_vol.placeVolume(c_vol, Transform3D(c_rot, c_pos));
} else if (x_rot) {
Position c_pos(0, 0, zoff);
pv = m_vol.placeVolume(c_vol, Transform3D(RotationZYX(x_rot.z(0), x_rot.y(0), x_rot.x(0)), c_pos));
} else if (x_pos) {
pv = m_vol.placeVolume(c_vol, Position(x_pos.x(0), x_pos.y(0), x_pos.z(0) + zoff));
} else {
pv = m_vol.placeVolume(c_vol, Position(0, 0, zoff));
}
c_vol.setRegion(description, x_comp.regionStr());
c_vol.setLimitSet(description, x_comp.limitsStr());
c_vol.setVisAttributes(description, x_comp.visStr());
if (x_comp.isSensitive()) {
pv.addPhysVolID("sensor", sensor_number++);
c_vol.setSensitiveDetector(sens);
sensitives[m_nam].push_back(pv);
module_thicknesses[m_nam] = {thickness_so_far + x_comp.thickness()/2.0, total_thickness-thickness_so_far - x_comp.thickness()/2.0};
// -------- create a measurement plane for the tracking surface attched to the sensitive volume -----
Vector3D u(-1., 0., 0.);
Vector3D v(0., -1., 0.);
Vector3D n(0., 0., 1.);
// Vector3D o( 0. , 0. , 0. ) ;
// compute the inner and outer thicknesses that need to be assigned to the tracking surface
// depending on wether the support is above or below the sensor
double inner_thickness = module_thicknesses[m_nam][0];
double outer_thickness = module_thicknesses[m_nam][1];
SurfaceType type(SurfaceType::Sensitive);
// if( isStripDetector )
// type.setProperty( SurfaceType::Measurement1D , true ) ;
VolPlane surf(c_vol, type, inner_thickness, outer_thickness, u, v, n); //,o ) ;
volplane_surfaces[m_nam].push_back(surf);
//--------------------------------------------
}
thickness_sum += x_comp.thickness();
thickness_so_far += x_comp.thickness();
// apply relative offsets in z-position used to stack components side-by-side
if (x_pos) {
thickness_sum += x_pos.z(0);
thickness_so_far += x_pos.z(0);
}
}
}
// now build the layers
for (xml_coll_t li(x_det, _U(layer)); li; ++li) {
xml_comp_t x_layer = li;
xml_comp_t x_barrel = x_layer.child(_U(barrel_envelope));
xml_comp_t x_layout = x_layer.child(_U(rphi_layout));
xml_comp_t z_layout = x_layer.child(_U(z_layout)); // Get the <z_layout> element.
int lay_id = x_layer.id();
string m_nam = x_layer.moduleStr();
string lay_nam = _toString(x_layer.id(), "layer%d");
Tube lay_tub(x_barrel.inner_r(), x_barrel.outer_r(), x_barrel.z_length() / 2.0);
Volume lay_vol(lay_nam, lay_tub, air); // Create the layer envelope volume.
lay_vol.setVisAttributes(description.visAttributes(x_layer.visStr()));
double phi0 = x_layout.phi0(); // Starting phi of first module.
double phi_tilt = x_layout.phi_tilt(); // Phi tilt of a module.
double rc = x_layout.rc(); // Radius of the module center.
int nphi = x_layout.nphi(); // Number of modules in phi.
double rphi_dr = x_layout.dr(); // The delta radius of every other module.
double phi_incr = (M_PI * 2) / nphi; // Phi increment for one module.
double phic = phi0; // Phi of the module center.
double z0 = z_layout.z0(); // Z position of first module in phi.
double nz = z_layout.nz(); // Number of modules to place in z.
double z_dr = z_layout.dr(); // Radial displacement parameter, of every other module.
Volume module_env = volumes[m_nam];
DetElement lay_elt(sdet, lay_nam, lay_id);
Placements& sensVols = sensitives[m_nam];
// the local coordinate systems of modules in dd4hep and acts differ
// see http://acts.web.cern.ch/ACTS/latest/doc/group__DD4hepPlugins.html
{
Acts::ActsExtension* layerExtension = new Acts::ActsExtension();
// layer is simple tube so no need to set envelope
layerExtension->addType("sensitive cylinder", "layer");
// Add the proto layer material
for(xml_coll_t lmat(x_layer, _Unicode(layer_material)); lmat; ++lmat) {
xml_comp_t x_layer_material = lmat;
xmlToProtoSurfaceMaterial(x_layer_material, *layerExtension, "layer_material");
}
lay_elt.addExtension<Acts::ActsExtension>(layerExtension);
}
// Z increment for module placement along Z axis.
// Adjust for z0 at center of module rather than
// the end of cylindrical envelope.
double z_incr = nz > 1 ? (2.0 * z0) / (nz - 1) : 0.0;
// Starting z for module placement along Z axis.
double module_z = -z0;
int module = 1;
// Loop over the number of modules in phi.
for (int ii = 0; ii < nphi; ii++) {
double dx = z_dr * std::cos(phic + phi_tilt); // Delta x of module position.
double dy = z_dr * std::sin(phic + phi_tilt); // Delta y of module position.
double x = rc * std::cos(phic); // Basic x module position.
double y = rc * std::sin(phic); // Basic y module position.
// Loop over the number of modules in z.
for (int j = 0; j < nz; j++) {
string module_name = _toString(module, "module%d");
DetElement mod_elt(lay_elt, module_name, module);
Transform3D tr(RotationZYX(0, ((M_PI / 2) - phic - phi_tilt), -M_PI / 2),
Position(x, y, module_z));
pv = lay_vol.placeVolume(module_env, tr);
pv.addPhysVolID("module", module);
mod_elt.setPlacement(pv);
for (size_t ic = 0; ic < sensVols.size(); ++ic) {
PlacedVolume sens_pv = sensVols[ic];
DetElement comp_de(mod_elt, std::string("de_") + sens_pv.volume().name(), module);
comp_de.setPlacement(sens_pv);
// ACTS extension
{
Acts::ActsExtension* sensorExtension = new Acts::ActsExtension();
//sensorExtension->addType("sensor", "detector");
comp_de.addExtension<Acts::ActsExtension>(sensorExtension);
}
//comp_de.setAttributes(description, sens_pv.volume(), x_layer.regionStr(), x_layer.limitsStr(),
// xml_det_t(xmleles[m_nam]).visStr());
//
volSurfaceList(comp_de)->push_back(volplane_surfaces[m_nam][ic]);
}
/// Increase counters etc.
module++;
// Adjust the x and y coordinates of the module.
x += dx;
y += dy;
// Flip sign of x and y adjustments.
dx *= -1;
dy *= -1;
// Add z increment to get next z placement pos.
module_z += z_incr;
}
phic += phi_incr; // Increment the phi placement of module.
rc += rphi_dr; // Increment the center radius according to dr parameter.
rphi_dr *= -1; // Flip sign of dr parameter.
module_z = -z0; // Reset the Z placement parameter for module.
}
// Create the PhysicalVolume for the layer.
pv = assembly.placeVolume(lay_vol); // Place layer in mother
pv.addPhysVolID("layer", lay_id); // Set the layer ID.
lay_elt.setAttributes(description, lay_vol, x_layer.regionStr(), x_layer.limitsStr(),
x_layer.visStr());
lay_elt.setPlacement(pv);
}
sdet.setAttributes(description, assembly, x_det.regionStr(), x_det.limitsStr(), x_det.visStr());
assembly.setVisAttributes(description.invisible());
pv = description.pickMotherVolume(sdet).placeVolume(assembly);
pv.addPhysVolID("system", det_id); // Set the subdetector system ID.
sdet.setPlacement(pv);
return sdet;
}
//@}
// clang-format off
DECLARE_DETELEMENT(BarrelTrackerWithFrame, create_BarrelTrackerWithFrame)
DECLARE_DETELEMENT(athena_TrackerBarrel, create_BarrelTrackerWithFrame)
DECLARE_DETELEMENT(athena_VertexBarrel, create_BarrelTrackerWithFrame)
DECLARE_DETELEMENT(athena_TOFBarrel, create_BarrelTrackerWithFrame)
src/CompositeTracker_geo.cpp 0 → 100644
View file @ 1008a184
//==========================================================================
// Based off DD4hep_SubDetectorAssembly
//
// This is a simple plugin to allow compositing different detectors
// into a single barrel or endcap for ACTS
// Note: positive/negative position strings differentiate between
// positive/negative endcaps
//--------------------------------------------------------------------------
#include "DD4hep/DetFactoryHelper.h"
#include "DD4hep/Printout.h"
#include "XML/Utilities.h"
#if defined(USE_ACTSDD4HEP)
#include "ActsDD4hep/ActsExtension.hpp"
#else
#include "Acts/Plugins/DD4hep/ActsExtension.hpp"
#endif
using namespace dd4hep;
using namespace dd4hep::detail;
static Ref_t create_element(Detector& description, xml_h e, Ref_t)
{
xml_det_t x_det(e);
const std::string det_name = x_det.nameStr();
DetElement sdet(det_name, x_det.id());
Volume vol;
Position pos;
const bool usePos = x_det.hasChild(_U(position));
sdet.setType("compound");
xml::setDetectorTypeFlag(e, sdet);
const std::string actsType = getAttrOrDefault(x_det, _Unicode(actsType), "endcap");
printout(DEBUG, det_name, "+++ Creating composite tracking detector (type: " + actsType + ")");
assert(actsType == "barrel" || actsType == "endcap");
// ACTS extension
{
Acts::ActsExtension* detWorldExt = new Acts::ActsExtension();
detWorldExt->addType(actsType, "detector");
sdet.addExtension<Acts::ActsExtension>(detWorldExt);
}
if (usePos) {
pos = Position(x_det.position().x(), x_det.position().y(), x_det.position().z());
}
vol = Assembly(det_name);
vol.setAttributes(description, x_det.regionStr(), x_det.limitsStr(), x_det.visStr());
Volume mother = description.pickMotherVolume(sdet);
PlacedVolume pv;
if (usePos) {
pv = mother.placeVolume(vol, pos);
} else {
pv = mother.placeVolume(vol);
}
sdet.setPlacement(pv);
for (xml_coll_t c(x_det, _U(composite)); c; ++c) {
xml_dim_t component = c;
const std::string nam = component.nameStr();
description.declareParent(nam, sdet);
}
return sdet;
}
DECLARE_DETELEMENT(athena_CompositeTracker, create_element)
src/CylinderTrackerBarrel_geo.cpp 0 → 100644
View file @ 1008a184
//==========================================================================
// Specialized generic detector constructor
//==========================================================================
#include "DD4hep/DetFactoryHelper.h"
#include "DD4hep/Printout.h"
#if defined(USE_ACTSDD4HEP)
#include "ActsDD4hep/ActsExtension.hpp"
#include "ActsDD4hep/ConvertMaterial.hpp"
#else
#include "Acts/Plugins/DD4hep/ActsExtension.hpp"
#include "Acts/Plugins/DD4hep/ConvertDD4hepMaterial.hpp"
#endif
using namespace std;
using namespace dd4hep;
using namespace dd4hep::detail;
/** A barrel tracker with a module that is curved (not flat).
*
* \ingroup tracking
*/
static Ref_t CylinderTrackerBarrel_create_detector(Detector& description, xml_h e, SensitiveDetector sens)
{
typedef vector<PlacedVolume> Placements;
xml_det_t x_det = e;
Material air = description.air();
int det_id = x_det.id();
string det_name = x_det.nameStr();
DetElement sdet(det_name, det_id);
//Acts::ActsExtension* barrelExtension = new Acts::ActsExtension();
//barrelExtension->addType("barrel", "detector");
//sdet.addExtension<Acts::ActsExtension>(barrelExtension);
Assembly assembly(det_name);
map<string, Volume> mod_volumes;
map<string, Placements> sensitives;
PlacedVolume pv;
sens.setType("tracker");
int n_modules = 0;
for (xml_coll_t mi(x_det, _U(module)); mi; ++mi) {
n_modules++;
xml_comp_t x_mod = mi;
xml_comp_t m_env = x_mod.child(_U(module_envelope));
string m_nam = x_mod.nameStr();
Assembly module_assembly(_toString(n_modules, "mod_assembly_%d"));
auto module_rmin = m_env.rmin();
auto module_thickness = m_env.thickness();
auto module_length = m_env.length();
auto module_phi = getAttrOrDefault(m_env, _Unicode(phi), 90.0);
Volume m_vol(
m_nam,
Tube(module_rmin, module_rmin + module_thickness, module_length / 2, -module_phi / 2.0, module_phi / 2.0), air);
int ncomponents = 0, sensor_number = 1;
module_assembly.placeVolume(m_vol, Position(-module_rmin, 0, 0));
mod_volumes[m_nam] = module_assembly;
m_vol.setVisAttributes(description.visAttributes(x_mod.visStr()));
auto comp_rmin = module_rmin;
for (xml_coll_t ci(x_mod, _U(module_component)); ci; ++ci, ++ncomponents) {
xml_comp_t x_comp = ci;
xml_comp_t x_pos = x_comp.position(false);
xml_comp_t x_rot = x_comp.rotation(false);
string c_nam = _toString(ncomponents, "component%d");
auto comp_thickness = x_comp.thickness();
comp_rmin = getAttrOrDefault(x_comp, _Unicode(rmin), comp_rmin);
auto comp_phi = getAttrOrDefault(x_comp, _Unicode(phi), module_phi);
auto comp_phi0 = getAttrOrDefault(x_comp, _Unicode(phi0), 0.0);
auto comp_length = getAttrOrDefault(x_comp, _Unicode(length), module_length);
Tube c_tube(comp_rmin, comp_rmin + comp_thickness, comp_length / 2, -comp_phi / 2.0 + comp_phi0,
comp_phi / 2.0 + comp_phi0);
Volume c_vol(c_nam, c_tube, description.material(x_comp.materialStr()));
PlacedVolume c_pv;
if (x_pos && x_rot) {
Position c_pos(x_pos.x(0), x_pos.y(0), x_pos.z(0));
RotationZYX c_rot(x_rot.z(0), x_rot.y(0), x_rot.x(0));
c_pv = m_vol.placeVolume(c_vol, Transform3D(c_rot, c_pos));
} else if (x_rot) {
c_pv = m_vol.placeVolume(c_vol, RotationZYX(x_rot.z(0), x_rot.y(0), x_rot.x(0)));
} else if (x_pos) {
c_pv = m_vol.placeVolume(c_vol, Position(x_pos.x(0), x_pos.y(0), x_pos.z(0)));
} else {
c_pv = m_vol.placeVolume(c_vol);
}
c_vol.setRegion(description, x_comp.regionStr());
c_vol.setLimitSet(description, x_comp.limitsStr());
c_vol.setVisAttributes(description, x_comp.visStr());
if (x_comp.isSensitive()) {
c_pv.addPhysVolID(_U(sensor), sensor_number++);
c_vol.setSensitiveDetector(sens);
sensitives[m_nam].push_back(c_pv);
}
comp_rmin = comp_rmin + comp_thickness;
}
}
for (xml_coll_t li(x_det, _U(layer)); li; ++li) {
xml_comp_t x_layer = li;
xml_comp_t x_barrel = x_layer.child(_U(barrel_envelope));
xml_comp_t x_layout = x_layer.child(_U(rphi_layout));
xml_comp_t z_layout = x_layer.child(_U(z_layout)); // Get the <z_layout> element.
int lay_id = x_layer.id();
string m_nam = x_layer.moduleStr();
string lay_nam = _toString(x_layer.id(), "layer%d");
Tube lay_tub(x_barrel.inner_r(), x_barrel.outer_r(), x_barrel.z_length() / 2);
Volume lay_vol(lay_nam, lay_tub, air); // Create the layer envelope volume.
lay_vol.setVisAttributes(description.visAttributes(x_layer.visStr()));
double phi0 = x_layout.phi0(); // Starting phi of first module.
double phi_tilt = x_layout.phi_tilt(); // Phi tilt of a module.
double rc = x_layout.rc(); // Radius of the module center.
int nphi = x_layout.nphi(); // Number of modules in phi.
double rphi_dr = x_layout.dr(); // The delta radius of every other module.
double phi_incr = (M_PI * 2) / nphi; // Phi increment for one module.
double phic = phi0; // Phi of the module center.
double z0 = z_layout.z0(); // Z position of first module in phi.
double nz = z_layout.nz(); // Number of modules to place in z.
double z_dr = z_layout.dr(); // Radial displacement parameter, of every other module.
Volume m_env = mod_volumes[m_nam];
DetElement lay_elt(sdet, _toString(x_layer.id(), "layer%d"), lay_id);
///Acts::ActsExtension* layerExtension = new Acts::ActsExtension();
///layerExtension->addType("sensitive cylinder", "layer");
/////// layerExtension->addValue(10. * Acts::UnitConstants::mm, "r", "envelope");
///lay_elt.addExtension<Acts::ActsExtension>(layerExtension);
Placements& sensVols = sensitives[m_nam];
// Z increment for module placement along Z axis.
// Adjust for z0 at center of module rather than
// the end of cylindrical envelope.
double z_incr = nz > 1 ? (2.0 * z0) / (nz - 1) : 0.0;
// Starting z for module placement along Z axis.
double module_z = -z0;
int module = 1;
// Loop over the number of modules in phi.
for (int ii = 0; ii < nphi; ii++) {
double dx = z_dr * std::cos(phic + phi_tilt); // Delta x of module position.
double dy = z_dr * std::sin(phic + phi_tilt); // Delta y of module position.
double x = rc * std::cos(phic); // Basic x module position.
double y = rc * std::sin(phic); // Basic y module position.
// Loop over the number of modules in z.
for (int j = 0; j < nz; j++) {
string module_name = _toString(module, "module%d");
DetElement mod_elt(lay_elt, module_name, module);
// Module PhysicalVolume.
// Transform3D
// tr(RotationZYX(0,-((M_PI/2)-phic-phi_tilt),M_PI/2),Position(x,y,module_z));
// NOTE (Nikiforos, 26/08 Rotations needed to be fixed so that component1 (silicon) is on the
// outside
Transform3D tr(RotationZYX(phic - phi_tilt, 0, 0), Position(x, y, module_z));
pv = lay_vol.placeVolume(m_env, tr);
pv.addPhysVolID("module", module);
mod_elt.setPlacement(pv);
for (size_t ic = 0; ic < sensVols.size(); ++ic) {
PlacedVolume sens_pv = sensVols[ic];
DetElement comp_elt(mod_elt, sens_pv.volume().name(), module);
comp_elt.setPlacement(sens_pv);
//Acts::ActsExtension* moduleExtension = new Acts::ActsExtension();
//comp_elt.addExtension<Acts::ActsExtension>(moduleExtension);
}
/// Increase counters etc.
module++;
// Adjust the x and y coordinates of the module.
x += dx;
y += dy;
// Flip sign of x and y adjustments.
dx *= -1;
dy *= -1;
// Add z increment to get next z placement pos.
module_z += z_incr;
}
phic += phi_incr; // Increment the phi placement of module.
rc += rphi_dr; // Increment the center radius according to dr parameter.
rphi_dr *= -1; // Flip sign of dr parameter.
module_z = -z0; // Reset the Z placement parameter for module.
}
// Create the PhysicalVolume for the layer.
pv = assembly.placeVolume(lay_vol); // Place layer in mother
pv.addPhysVolID("layer", lay_id); // Set the layer ID.
lay_elt.setAttributes(description, lay_vol, x_layer.regionStr(), x_layer.limitsStr(), x_layer.visStr());
lay_elt.setPlacement(pv);
}
sdet.setAttributes(description, assembly, x_det.regionStr(), x_det.limitsStr(), x_det.visStr());
// assembly.setVisAttributes(description.invisible());
pv = description.pickMotherVolume(sdet).placeVolume(assembly);
pv.addPhysVolID("system", det_id); // Set the subdetector system ID.
pv.addPhysVolID("barrel", 1); // Flag this as a barrel subdetector.
sdet.setPlacement(pv);
return sdet;
}
// clang-format off
DECLARE_DETELEMENT(athena_CylinderTrackerBarrel, CylinderTrackerBarrel_create_detector)
DECLARE_DETELEMENT(athena_MMTrackerBarrel, CylinderTrackerBarrel_create_detector)
DECLARE_DETELEMENT(athena_RWellTrackerBarrel, CylinderTrackerBarrel_create_detector)
DECLARE_DETELEMENT(athena_CylinderVertexBarrel, CylinderTrackerBarrel_create_detector)
src/DIRC_geo.cpp 0 → 100644
View file @ 1008a184
#include "DD4hep/DetFactoryHelper.h"
#include "DD4hep/OpticalSurfaces.h"
#include "DD4hep/Printout.h"
#include "DDRec/DetectorData.h"
#include "DDRec/Surface.h"
#include <XML/Helper.h>
//////////////////////////////////
// Central Barrel DIRC
//////////////////////////////////
using namespace std;
using namespace dd4hep;
// Fixed Trap constructor. This function is a workaround of this bug:
// https://github.com/AIDASoft/DD4hep/issues/850
// Should be used instead of dd4hep::Trap(pName, pZ, pY, pX, pLTX) constructor
dd4hep::Trap MakeTrap( const std::string& pName, double pZ, double pY, double pX, double pLTX );
static Ref_t createDetector(Detector& desc, xml_h e, SensitiveDetector sens)
{
xml_det_t xml_det = e;
string det_name = xml_det.nameStr();
int det_id = xml_det.id();
// Main detector xml element
xml_dim_t dirc_dim = xml_det.dimensions();
xml_dim_t dirc_pos = xml_det.position();
xml_dim_t dirc_rot = xml_det.rotation();
double det_rin = dirc_dim.rmin();
double det_rout = dirc_dim.rmax();
double SizeZ = dirc_dim.length();
// DEBUG
// double mirror_r1 = x_det.attr<double>(_Unicode(r1));
// DIRC box:
xml_comp_t xml_box_module = xml_det.child(_U(module));
Material Vacuum = desc.material("Vacuum");
Material air = desc.material("AirOptical");
Material quartz = desc.material("Quartz");
Material epotek = desc.material("Epotek");
Material nlak33a = desc.material("Nlak33a");
auto& bar_material = quartz;
auto mirror_material = desc.material("Aluminum"); // mirror material
Tube det_geo(det_rin, det_rout, SizeZ / 2., 0., 360.0 * deg);
//Volume det_volume("DIRC", det_geo, Vacuum);
Assembly det_volume("DIRC");
det_volume.setVisAttributes(desc.visAttributes(xml_det.visStr()));
DetElement det(det_name, det_id);
Volume mother_vol = desc.pickMotherVolume(det);
Transform3D tr(RotationZYX(0, dirc_rot.theta(), 0.0), Position(0.0, 0.0, dirc_pos.z()));
PlacedVolume det_plvol = mother_vol.placeVolume(det_volume, tr);
det_plvol.addPhysVolID("system", det_id);
det.setPlacement(det_plvol);
// Parts Dimentions
int fLensId = 6; // focusing system
// 0 no lens
// 1 spherical lens
// 3 3-layer spherical lens
// 6 3-layer cylindrical lens
// 10 ideal lens (thickness = 0, ideal focusing)
int fGeomType = 0; // Full DIRC - 0, 1 only one plate
int fRunType = 0; // 0, 10 - simulation, 1, 5 - lookup table, 2,3,4 - reconstruction
double fPrizm[4];
fPrizm[0] = 360 * mm;
fPrizm[1] = 300 * mm;
fPrizm[3] = 50 * mm;
fPrizm[2] = fPrizm[3] + 300 * tan(32 * deg) * mm;
double fBarsGap = 0.15 * mm;
std::cout << "DIRC: fPrizm[2] " << fPrizm[2] << std::endl;
double fdTilt = 80 * deg;
double fPrizmT[6];
fPrizmT[0] = 390 * mm;
fPrizmT[1] = (400 - 290 * cos(fdTilt)) * mm; //
fPrizmT[2] = 290 * sin(fdTilt) * mm; // hight
fPrizmT[3] = 50 * mm; // face
fPrizmT[4] = 290 * mm;
fPrizmT[5] = 290 * cos(fdTilt)* mm;
double fMirror[3];
fMirror[0] = 20 * mm;
fMirror[1] = fPrizm[0];
fMirror[2] = 1 * mm;
// fPrizm[0] = 170; fPrizm[1] = 300; fPrizm[2] = 50+300*tan(45*deg); fPrizm[3] = 50;
// double fBar[3];
// fBar[0] = 17 * mm;
// fBar[1] = (fPrizm[0] - (fNBar - 1) * fBarsGap) / fNBar;
// fBar[2] = 1050 * mm; // 4200; //4200
double fMcpTotal[3];
double fMcpActive[3];
fMcpTotal[0] = fMcpTotal[1] = 53 + 4;
fMcpTotal[2] = 1*mm;
fMcpActive[0] = fMcpActive[1] = 53;
fMcpActive[2] = 1*mm;
double fLens[4];
fLens[0] = fLens[1] = 40 * mm;
fLens[2] = 10 * mm;
double fRadius = (det_rin + det_rout)/2;
double fBoxWidth = fPrizm[0];
double fFd[3];
fFd[0] = fBoxWidth;
fFd[1] = fPrizm[2];
fFd[2] = 1 * mm;
fLens[0] = fPrizm[3];
fLens[1] = fPrizm[0];
fLens[2] = 12 * mm;
// Getting box XML
const int fNBoxes = xml_box_module.repeat();
const double box_width = xml_box_module.width();
const double box_height = xml_box_module.height();
const double box_length = xml_box_module.length() + 550*mm;
// The DIRC
Assembly dirc_module("DIRCModule");
//Volume lDirc("lDirc", gDirc, air);
dirc_module.setVisAttributes(desc.visAttributes(xml_box_module.visStr()));
// FD... whatever F and D is
xml_comp_t xml_fd = xml_box_module.child(_Unicode(fd));
Box gFd("gFd", xml_fd.height()/2, xml_fd.width()/2, xml_fd.thickness()/2);
Volume lFd ("lFd", gFd, desc.material(xml_fd.materialStr()));
lFd.setVisAttributes(desc.visAttributes(xml_fd.visStr()));
//lFd.setSensitiveDetector(sens);
// The Bar
xml_comp_t xml_bar = xml_box_module.child(_Unicode(bar));
double bar_height = xml_bar.height();
double bar_width = xml_bar.width();
double bar_length = xml_bar.length();
Box gBar("gBar", bar_height/2, bar_width/2, bar_length/2);
Volume lBar("lBar", gBar, desc.material(xml_bar.materialStr()));
lBar.setVisAttributes(desc.visAttributes(xml_bar.visStr()));
// Glue
xml_comp_t xml_glue = xml_box_module.child(_Unicode(glue));
double glue_thickness = xml_glue.thickness(); // 0.05 * mm;
Box gGlue("gGlue", bar_height/2, bar_width/2, glue_thickness/2);
Volume lGlue("lGlue", gGlue, desc.material(xml_glue.materialStr()));
lGlue.setVisAttributes(desc.visAttributes(xml_glue.visStr()));
sens.setType("tracker");
lBar.setSensitiveDetector(sens);
int bars_repeat_z = 4; // TODO parametrize!
double bar_assm_length = (bar_length + glue_thickness) * bars_repeat_z;
int fNBar = xml_bar.repeat();
double bar_gap = xml_bar.gap();
for (int y_index = 0; y_index < fNBar; y_index++) {
double shift_y = y_index * (bar_width + bar_gap) - 0.5 * box_width + 0.5 * bar_width;
for (int z_index = 0; z_index < bars_repeat_z; z_index++) {
double z = -0.5 * bar_assm_length + 0.5 * bar_length + (bar_length + glue_thickness) * z_index;
auto placed_bar = dirc_module.placeVolume(lBar, Position(0, shift_y, z));
dirc_module.placeVolume(lGlue, Position(0, shift_y, z + 0.5 * (bar_length + glue_thickness)));
placed_bar.addPhysVolID("section", z_index);
placed_bar.addPhysVolID("bar", y_index);
}
}
// The Mirror
xml_comp_t xml_mirror = xml_box_module.child(_Unicode(mirror));
Box gMirror("gMirror", xml_mirror.height()/2, xml_mirror.width()/2, xml_mirror.thickness()/2);
Volume lMirror("lMirror", gMirror, desc.material(xml_mirror.materialStr()));
dirc_module.placeVolume(lMirror, Position(0, 0, -0.5 * (bar_assm_length - xml_mirror.thickness())));
lMirror.setVisAttributes(desc.visAttributes(xml_mirror.visStr()));
// The mirror optical surface
OpticalSurfaceManager surfMgr = desc.surfaceManager();
auto surf = surfMgr.opticalSurface("MirrorOpticalSurface");
SkinSurface skin(desc, det, Form("dirc_mirror_optical_surface"), surf, lMirror);
skin.isValid();
// LENS
// Lens volumes
Volume lLens1;
Volume lLens2;
Volume lLens3;
double lensMinThikness = 2.0 * mm;
double layer12 = lensMinThikness * 2;
// r1 = (r1==0)? 27.45: r1;
// r2 = (r2==0)? 20.02: r2;
double r1 = 33 * mm;
double r2 = 24 * mm;
double shight = 25 * mm;
Position zTrans1(0, 0, -r1 - fLens[2] / 2. + r1 - sqrt(r1 * r1 - shight / 2. * shight / 2.) + lensMinThikness);
Position zTrans2(0, 0, -r2 - fLens[2] / 2. + r2 - sqrt(r2 * r2 - shight / 2. * shight / 2.) + layer12);
Box gfbox("fbox", 0.5 * fLens[0], 0.5 * fLens[1], 0.5 * fLens[2]);
Box gcbox("cbox", 0.5 * fLens[0], 0.5 * fLens[1] + 1*mm, 0.5 * fLens[2]);
// Volume gfbox_volume("gfbox_volume", gfbox, bar_material);
// lDirc.placeVolume(gfbox_volume, Position(0, 0, 0));
//
// Volume gcbox_volume("gcbox_volume", gcbox, bar_material);
// lDirc.placeVolume(gcbox_volume, Position(0, 0, 50));
Position tTrans1(0.5 * (fLens[0] + shight), 0, -fLens[2] + layer12);
Position tTrans0(-0.5 * (fLens[0] + shight), 0, -fLens[2] + layer12);
SubtractionSolid tubox("tubox", gfbox, gcbox, tTrans1);
SubtractionSolid gubox("gubox", tubox, gcbox, tTrans0);
// Volume tubox_volume("tubox_volume", tubox, bar_material);
// lDirc.placeVolume(tubox_volume, Position(0, 0, 100));
//
// Volume gubox_volume("gubox_volume", gubox, bar_material);
// lDirc.placeVolume(gubox_volume, Position(0, 0, 150));
Tube gcylinder1("Cylinder1", 0, r1, 0.5 * fLens[1], 0 * deg, 360 * deg);
Tube gcylinder2("Cylinder2", 0, r2, 0.5 * fLens[1] - 0.5*mm, 0 * deg, 360 * deg);
Tube gcylinder1c("Cylinder1c", 0, r1, 0.5 * fLens[1] + 0.5*mm, 0 * deg, 360 * deg);
Tube gcylinder2c("Cylinder2c", 0, r2, 0.5 * fLens[1] + 0.5*mm, 0 * deg, 360 * deg);
RotationX xRot(-M_PI / 2.);
IntersectionSolid gLens1("Lens1", gubox, gcylinder1, Transform3D(xRot, zTrans1));
SubtractionSolid gLenst("temp", gubox, gcylinder1c, Transform3D(xRot, zTrans1));
// Volume gLens1_volume("gLens1_volume", gLens1, bar_material);
// lDirc.placeVolume(gLens1_volume, Position(0, 0, 200));
//
// Volume gLenst_volume("gLenst_volume", gLenst, bar_material);
// lDirc.placeVolume(gLenst_volume, Position(0, 0, 250));
IntersectionSolid gLens2("Lens2", gLenst, gcylinder2, Transform3D(xRot, zTrans2));
SubtractionSolid gLens3("Lens3", gLenst, gcylinder2c, Transform3D(xRot, zTrans2));
lLens1 = Volume("lLens1", gLens1, bar_material);
lLens2 = Volume("lLens2", gLens2, nlak33a);
lLens3 = Volume("lLens3", gLens3, bar_material);
lLens1.setVisAttributes(desc.visAttributes("DIRCLens1"));
lLens2.setVisAttributes(desc.visAttributes("DIRCLens2"));
lLens3.setVisAttributes(desc.visAttributes("DIRCLens3"));
double shifth = 0.5 * (bar_assm_length + fLens[2]);
// fmt::print("LENS HERE shifth={}\n", shifth);
lLens1.setVisAttributes(desc.visAttributes("AnlTeal"));
dirc_module.placeVolume(lLens1, Position(0, 0, shifth));
dirc_module.placeVolume(lLens2, Position(0, 0, shifth));
dirc_module.placeVolume(lLens3, Position(0, 0, shifth));
// The Prizm
Trap gPrizm = MakeTrap("gPrizm", fPrizm[0], fPrizm[1], fPrizm[2], fPrizm[3]);
Volume lPrizm("lPrizm", gPrizm, bar_material);
lPrizm.setVisAttributes(desc.visAttributes("DIRCPrism"));
//G4RotationMatrix *fdRot = new G4RotationMatrix();
//G4RotationMatrix *fdrot = new G4RotationMatrix();
double evshiftz = 0.5 * bar_assm_length + fPrizm[1] + fMcpActive[2] / 2. + fLens[2];
double evshiftx = -3*mm;
double prism_shift_x = (fPrizm[2] + fPrizm[3]) / 4. - 0.5 * fPrizm[3] + 1.5*mm;
double prism_shift_z = 0.5 * (bar_assm_length + fPrizm[1]) + fLens[2];
Position fPrismShift(prism_shift_x, 0, prism_shift_z);
dirc_module.placeVolume(lPrizm, Transform3D(xRot, fPrismShift));
dirc_module.placeVolume(lFd, Position(0.5 * fFd[1] - 0.5 * fPrizm[3] - evshiftx, 0, evshiftz));
double dphi = 2 * M_PI / (double)fNBoxes;
for (int i = 0; i < fNBoxes; i++) {
double phi = dphi * i;
double dx = -fRadius * cos(phi);
double dy = -fRadius * sin(phi);
//G4RotationMatrix *tRot = new G4RotationMatrix();
Transform3D tr(RotationZ(phi+M_PI), Position(dx, dy, 0));
PlacedVolume box_placement = det_volume.placeVolume(dirc_module, tr);
box_placement.addPhysVolID("module", i);
// fmt::print("placing dircbox # {} -tphi={:.0f} dx={:.0f}, dy={:.0f}\n", i, phi/deg, dx/cm, dy/cm);
//new G4PVPlacement(tRot, G4ThreeVector(dx, dy, 0), lDirc, "wDirc", lExpHall, false, i);
}
//////////////////
// DIRC Bars
//////////////////
// double bar_radius = 83.65 * cm;
// double bar_length = SizeZ;
// double bar_width = 42. * cm;
// double bar_thicknes = 1.7 * cm;
// int bar_count = 2 * M_PI * bar_radius / bar_width;
// double bar_dphi = 2 * 3.1415926 / bar_count;
// Material bar_material = desc.material("Quartz");
// Box bar_geo(bar_thicknes / 2., bar_width / 2., bar_length / 2.);
// Volume bar_volume("cb_DIRC_bars_Logix", bar_geo, bar_material);
// bar_volume.setVisAttributes(desc.visAttributes(xml_det.visStr()));
// sens.setType("tracker");
// bar_volume.setSensitiveDetector(sens);
// for (int ia = 0; ia < bar_count; ia++) {
// double phi = (ia * (bar_dphi));
// double x = -bar_radius * cos(phi);
// double y = -bar_radius * sin(phi);
// Transform3D tr(RotationZ(bar_dphi * ia), Position(x, y, 0));
// PlacedVolume barPV = det_volume.placeVolume(bar_volume, tr);
// barPV.addPhysVolID("module", ia);
// }
return det;
}
DECLARE_DETELEMENT(cb_DIRC, createDetector)
dd4hep::Trap MakeTrap( const std::string& pName, double pZ, double pY, double pX, double pLTX )
{
// Fixed Trap constructor. This function is a workaround of this bug:
// https://github.com/AIDASoft/DD4hep/issues/850
// Should be used instead of dd4hep::Trap(pName, pZ, pY, pX, pLTX) constructor
double fDz = 0.5*pZ;
double fTthetaCphi = 0;
double fTthetaSphi = 0;
double fDy1 = 0.5*pY;
double fDx1 = 0.5*pX;
double fDx2 = 0.5*pLTX;
double fTalpha1 = 0.5*(pLTX - pX)/pY;
double fDy2 = fDy1;
double fDx3 = fDx1;
double fDx4 = fDx2;
double fTalpha2 = fTalpha1;
return Trap(pName, fDz, fTthetaCphi, fTthetaSphi,
fDy1, fDx1, fDx2, fTalpha1,
fDy2, fDx3, fDx4, fTalpha2);
}
src/DRICH_geo.cpp 0 → 100644
View file @ 1008a184
//==========================================================================
// dRICh: Dual Ring Imaging Cherenkov Detector
//--------------------------------------------------------------------------
//
// Author: Christopher Dilks (Duke University)
//
// - Design Adapted from Standalone Fun4all and GEMC implementations
// [ Evaristo Cisbani, Cristiano Fanelli, Alessio Del Dotto, et al. ]
//
//==========================================================================
#include <XML/Helper.h>
#include "TMath.h"
#include "TString.h"
#include "GeometryHelpers.h"
#include "Math/Point2D.h"
#include "DDRec/Surface.h"
#include "DDRec/DetectorData.h"
#include "DD4hep/OpticalSurfaces.h"
#include "DD4hep/DetFactoryHelper.h"
#include "DD4hep/Printout.h"
using namespace dd4hep;
using namespace dd4hep::rec;
// create the detector
static Ref_t createDetector(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);
xml::Component dims = detElem.dimensions();
OpticalSurfaceManager surfMgr = desc.surfaceManager();
// attributes -----------------------------------------------------------
// - vessel
double vesselZmin = dims.attr<double>(_Unicode(zmin));
double vesselLength = dims.attr<double>(_Unicode(length));
double vesselRmin0 = dims.attr<double>(_Unicode(rmin0));
double vesselRmin1 = dims.attr<double>(_Unicode(rmin1));
double vesselRmax0 = dims.attr<double>(_Unicode(rmax0));
double vesselRmax1 = dims.attr<double>(_Unicode(rmax1));
double vesselRmax2 = dims.attr<double>(_Unicode(rmax2));
double snoutLength = dims.attr<double>(_Unicode(snout_length));
int nSectors = dims.attr<int>(_Unicode(nsectors));
double wallThickness = dims.attr<double>(_Unicode(wall_thickness));
double windowThickness = dims.attr<double>(_Unicode(window_thickness));
auto vesselMat = desc.material(detElem.attr<std::string>(_Unicode(material)));
auto gasvolMat = desc.material(detElem.attr<std::string>(_Unicode(gas)));
auto vesselVis = desc.visAttributes(detElem.attr<std::string>(_Unicode(vis_vessel)));
auto gasvolVis = desc.visAttributes(detElem.attr<std::string>(_Unicode(vis_gas)));
// - radiator (applies to aerogel and filter)
auto radiatorElem = detElem.child(_Unicode(radiator));
double radiatorRmin = radiatorElem.attr<double>(_Unicode(rmin));
double radiatorRmax = radiatorElem.attr<double>(_Unicode(rmax));
double radiatorPhiw = radiatorElem.attr<double>(_Unicode(phiw));
double radiatorPitch = radiatorElem.attr<double>(_Unicode(pitch));
double radiatorFrontplane = radiatorElem.attr<double>(_Unicode(frontplane));
// - aerogel
auto aerogelElem = radiatorElem.child(_Unicode(aerogel));
auto aerogelMat = desc.material(aerogelElem.attr<std::string>(_Unicode(material)));
auto aerogelVis = desc.visAttributes(aerogelElem.attr<std::string>(_Unicode(vis)));
double aerogelThickness = aerogelElem.attr<double>(_Unicode(thickness));
// - filter
auto filterElem = radiatorElem.child(_Unicode(filter));
auto filterMat = desc.material(filterElem.attr<std::string>(_Unicode(material)));
auto filterVis = desc.visAttributes(filterElem.attr<std::string>(_Unicode(vis)));
double filterThickness = filterElem.attr<double>(_Unicode(thickness));
// - mirror
auto mirrorElem = detElem.child(_Unicode(mirror));
auto mirrorMat = desc.material(mirrorElem.attr<std::string>(_Unicode(material)));
auto mirrorVis = desc.visAttributes(mirrorElem.attr<std::string>(_Unicode(vis)));
auto mirrorSurf = surfMgr.opticalSurface(mirrorElem.attr<std::string>(_Unicode(surface)));
double mirrorBackplane = mirrorElem.attr<double>(_Unicode(backplane));
double mirrorThickness = mirrorElem.attr<double>(_Unicode(thickness));
double mirrorRmin = mirrorElem.attr<double>(_Unicode(rmin));
double mirrorRmax = mirrorElem.attr<double>(_Unicode(rmax));
double mirrorPhiw = mirrorElem.attr<double>(_Unicode(phiw));
double focusTuneZ = mirrorElem.attr<double>(_Unicode(focus_tune_z));
double focusTuneX = mirrorElem.attr<double>(_Unicode(focus_tune_x));
// - sensor module
auto sensorElem = detElem.child(_Unicode(sensors)).child(_Unicode(module));
auto sensorMat = desc.material(sensorElem.attr<std::string>(_Unicode(material)));
auto sensorVis = desc.visAttributes(sensorElem.attr<std::string>(_Unicode(vis)));
auto sensorSurf = surfMgr.opticalSurface(sensorElem.attr<std::string>(_Unicode(surface)));
double sensorSide = sensorElem.attr<double>(_Unicode(side));
double sensorGap = sensorElem.attr<double>(_Unicode(gap));
double sensorThickness = sensorElem.attr<double>(_Unicode(thickness));
// - sensor sphere
auto sensorSphElem = detElem.child(_Unicode(sensors)).child(_Unicode(sphere));
double sensorSphRadius = sensorSphElem.attr<double>(_Unicode(radius));
double sensorSphCenterX = sensorSphElem.attr<double>(_Unicode(centerx));
double sensorSphCenterZ = sensorSphElem.attr<double>(_Unicode(centerz));
// - sensor sphere patch cuts
auto sensorSphPatchElem = detElem.child(_Unicode(sensors)).child(_Unicode(sphericalpatch));
double sensorSphPatchPhiw = sensorSphPatchElem.attr<double>(_Unicode(phiw));
double sensorSphPatchRmin = sensorSphPatchElem.attr<double>(_Unicode(rmin));
double sensorSphPatchRmax = sensorSphPatchElem.attr<double>(_Unicode(rmax));
double sensorSphPatchZmin = sensorSphPatchElem.attr<double>(_Unicode(zmin));
// - debugging switches
int debug_optics_mode = detElem.attr<int>(_Unicode(debug_optics));
bool debug_mirror = mirrorElem.attr<bool>(_Unicode(debug));
bool debug_sensors = sensorSphElem.attr<bool>(_Unicode(debug));
// if debugging optics, override some settings
bool debug_optics = debug_optics_mode > 0;
if(debug_optics) {
printout(WARNING,"DRich_geo","DEBUGGING DRICH OPTICS");
switch(debug_optics_mode) {
case 1: vesselMat = aerogelMat = filterMat = sensorMat = gasvolMat = desc.material("VacuumOptical"); break;
case 2: vesselMat = aerogelMat = filterMat = sensorMat = desc.material("VacuumOptical"); break;
default: printout(FATAL,"DRich_geo","UNKNOWN debug_optics_mode"); return det;
};
aerogelVis = sensorVis = mirrorVis;
gasvolVis = vesselVis = desc.invisible();
};
// BUILD VESSEL ====================================================================
/* - `vessel`: aluminum enclosure, the mother volume of the dRICh
* - `gasvol`: gas volume, which fills `vessel`; all other volumes defined below
* are children of `gasvol`
* - the dRICh vessel geometry has two regions: the snout refers to the conic region
* in the front, housing the aerogel, while the tank refers to the cylindrical
* region, housing the rest of the detector components
*/
// derived attributes
double tankLength = vesselLength - snoutLength;
double vesselZmax = vesselZmin + vesselLength;
// snout solids
double boreDelta = vesselRmin1 - vesselRmin0;
double snoutDelta = vesselRmax1 - vesselRmax0;
Cone vesselSnout(
snoutLength/2.0,
vesselRmin0,
vesselRmax0,
vesselRmin0 + boreDelta * snoutLength / vesselLength,
vesselRmax1
);
Cone gasvolSnout(
/* note: `gasvolSnout` extends a bit into the tank, so it touches `gasvolTank`
* - the extension distance is equal to the tank `windowThickness`, so the
* length of `gasvolSnout` == length of `vesselSnout`
* - the extension backplane radius is calculated using similar triangles
*/
snoutLength/2.0,
vesselRmin0 + wallThickness,
vesselRmax0 - wallThickness,
vesselRmin0 + boreDelta * (snoutLength-windowThickness) / vesselLength + wallThickness,
vesselRmax1 - wallThickness + windowThickness * (vesselRmax1 - vesselRmax0) / snoutLength
);
// tank solids
Cone vesselTank(
tankLength/2.0,
vesselSnout.rMin2(),
vesselRmax2,
vesselRmin1,
vesselRmax2
);
Cone gasvolTank(
tankLength/2.0 - windowThickness,
gasvolSnout.rMin2(),
vesselRmax2 - wallThickness,
vesselRmin1 + wallThickness,
vesselRmax2 - wallThickness
);
// snout + tank solids
UnionSolid vesselUnion(
vesselTank,
vesselSnout,
Position(0., 0., -vesselLength/2.)
);
UnionSolid gasvolUnion(
gasvolTank,
gasvolSnout,
Position(0., 0., -vesselLength/2. + windowThickness)
);
// extra solids for `debug_optics` only
Box vesselBox(1001,1001,1001);
Box gasvolBox(1000,1000,1000);
// choose vessel and gasvol solids (depending on `debug_optics_mode` (0=disabled))
Solid vesselSolid, gasvolSolid;
switch(debug_optics_mode) {
case 0: vesselSolid=vesselUnion; gasvolSolid=gasvolUnion; break; // `!debug_optics`
case 1: vesselSolid=vesselBox; gasvolSolid=gasvolBox; break;
case 2: vesselSolid=vesselBox; gasvolSolid=gasvolUnion; break;
};
// volumes
Volume vesselVol(detName, vesselSolid, vesselMat);
Volume gasvolVol(detName+"_gas", gasvolSolid, gasvolMat);
vesselVol.setVisAttributes(vesselVis);
gasvolVol.setVisAttributes(gasvolVis);
// reference positions
// - the vessel is created such that the center of the cylindrical tank volume
// coincides with the origin; this is called the "origin position" of the vessel
// - when the vessel (and its children volumes) is placed, it is translated in
// the z-direction to be in the proper ATHENA-integration location
// - these reference positions are for the frontplane and backplane of the vessel,
// with respect to the vessel origin position
auto originFront = Position(0., 0., -tankLength/2.0 - snoutLength );
auto originBack = Position(0., 0., tankLength/2.0 );
// initialize sensor centroids (used for mirror parameterization below); this is
// the average (x,y,z) of the placed sensors, w.r.t. originFront
double sensorCentroidX = 0;
double sensorCentroidZ = 0;
int sensorCount = 0;
// sensitive detector type
sens.setType("tracker");
// BUILD RADIATOR ====================================================================
// attributes
double airGap = 0.01*mm; // air gap between aerogel and filter (FIXME? actually it's currently a gas gap)
// solid and volume: create aerogel and filter
Cone aerogelSolid(
aerogelThickness/2,
radiatorRmin,
radiatorRmax,
radiatorRmin + boreDelta * aerogelThickness / vesselLength,
radiatorRmax + snoutDelta * aerogelThickness / snoutLength
);
Cone filterSolid(
filterThickness/2,
radiatorRmin + boreDelta * (aerogelThickness + airGap) / vesselLength,
radiatorRmax + snoutDelta * (aerogelThickness + airGap) / snoutLength,
radiatorRmin + boreDelta * (aerogelThickness + airGap + filterThickness) / vesselLength,
radiatorRmax + snoutDelta * (aerogelThickness + airGap + filterThickness) / snoutLength
);
Volume aerogelVol( detName+"_aerogel", aerogelSolid, aerogelMat );
Volume filterVol( detName+"_filter", filterSolid, filterMat );
aerogelVol.setVisAttributes(aerogelVis);
filterVol.setVisAttributes(filterVis);
// aerogel placement and surface properties
// TODO [low-priority]: define skin properties for aerogel and filter
auto radiatorPos = Position(0., 0., radiatorFrontplane) + originFront;
auto aerogelPV = gasvolVol.placeVolume(aerogelVol,
Translation3D(radiatorPos.x(), radiatorPos.y(), radiatorPos.z()) // re-center to originFront
* RotationY(radiatorPitch) // change polar angle to specified pitch // FIXME: probably broken (not yet in use anyway)
);
DetElement aerogelDE(det, "aerogel_de", 0);
aerogelDE.setPlacement(aerogelPV);
//SkinSurface aerogelSkin(desc, aerogelDE, "mirror_optical_surface", aerogelSurf, aerogelVol);
//aerogelSkin.isValid();
// filter placement and surface properties
if(!debug_optics) {
auto filterPV = gasvolVol.placeVolume(filterVol,
Translation3D(0., 0., airGap) // add an air gap
* Translation3D(radiatorPos.x(), radiatorPos.y(), radiatorPos.z()) // re-center to originFront
* RotationY(radiatorPitch) // change polar angle
* Translation3D(0., 0., (aerogelThickness+filterThickness)/2.) // move to aerogel backplane
);
DetElement filterDE(det, "filter_de", 0);
filterDE.setPlacement(filterPV);
//SkinSurface filterSkin(desc, filterDE, "mirror_optical_surface", filterSurf, filterVol);
//filterSkin.isValid();
};
// SECTOR LOOP //////////////////////////////////
for(int isec=0; isec<nSectors; isec++) {
// debugging filters, limiting the number of sectors
if( (debug_mirror||debug_sensors||debug_optics) && isec!=0) continue;
// sector rotation about z axis
double sectorRotation = isec * 360/nSectors * degree;
std::string secName = "sec" + std::to_string(isec);
// BUILD SENSORS ====================================================================
// if debugging sphere properties, restrict number of sensors drawn
if(debug_sensors) { sensorSide = 2*M_PI*sensorSphRadius / 64; };
// solid and volume: single sensor module
Box sensorSolid(sensorSide/2., sensorSide/2., sensorThickness/2.);
Volume sensorVol(detName+"_sensor_"+secName, sensorSolid, sensorMat);
sensorVol.setVisAttributes(sensorVis);
auto sensorSphPos = Position(sensorSphCenterX, 0., sensorSphCenterZ) + originFront;
// sensitivity
if(!debug_optics) sensorVol.setSensitiveDetector(sens);
// SENSOR MODULE LOOP ------------------------
/* ALGORITHM: generate sphere of positions
* - NOTE: there are two coordinate systems here:
* - "global" the main ATHENA coordinate system
* - "generator" (vars end in `Gen`) is a local coordinate system for
* generating points on a sphere; it is related to the global system by
* a rotation; we do this so the "patch" (subset of generated
* positions) of sensors we choose to build is near the equator, where
* point distribution is more uniform
* - PROCEDURE: loop over `thetaGen`, with subloop over `phiGen`, each divided evenly
* - the number of points to generate depends how many sensors (+`sensorGap`)
* can fit within each ring of constant `thetaGen` or `phiGen`
* - we divide the relevant circumference by the sensor
* size(+`sensorGap`), and this number is allowed to be a fraction,
* because likely we don't care about generating a full sphere and
* don't mind a "seam" at the overlap point
* - if we pick a patch of the sphere near the equator, and not near
* the poles or seam, the sensor distribution will appear uniform
*/
// initialize module number for this sector
int imod=0;
// thetaGen loop: iterate less than "0.5 circumference / sensor size" times
double nTheta = M_PI*sensorSphRadius / (sensorSide+sensorGap);
for(int t=0; t<(int)(nTheta+0.5); t++) {
double thetaGen = t/((double)nTheta) * M_PI;
// phiGen loop: iterate less than "circumference at this latitude / sensor size" times
double nPhi = 2*M_PI * sensorSphRadius * std::sin(thetaGen) / (sensorSide+sensorGap);
for(int p=0; p<(int)(nPhi+0.5); p++) {
double phiGen = p/((double)nPhi) * 2*M_PI - M_PI; // shift to [-pi,pi]
// determine global phi and theta
// - convert {radius,thetaGen,phiGen} -> {xGen,yGen,zGen}
double xGen = sensorSphRadius * std::sin(thetaGen) * std::cos(phiGen);
double yGen = sensorSphRadius * std::sin(thetaGen) * std::sin(phiGen);
double zGen = sensorSphRadius * std::cos(thetaGen);
// - convert {xGen,yGen,zGen} -> global {x,y,z} via rotation
double x = zGen;
double y = xGen;
double z = yGen;
// - convert global {x,y,z} -> global {phi,theta}
double phi = std::atan2(y,x);
double theta = std::acos(z/sensorSphRadius);
// shift global coordinates so we can apply spherical patch cuts
double zCheck = z + sensorSphCenterZ;
double xCheck = x + sensorSphCenterX;
double yCheck = y;
double rCheck = std::hypot(xCheck,yCheck);
double phiCheck = std::atan2(yCheck,xCheck);
// patch cut
bool patchCut =
std::fabs(phiCheck) < sensorSphPatchPhiw
&& zCheck > sensorSphPatchZmin
&& rCheck > sensorSphPatchRmin
&& rCheck < sensorSphPatchRmax;
if(debug_sensors) patchCut = std::fabs(phiCheck) < sensorSphPatchPhiw;
if(patchCut) {
// append sensor position to centroid calculation
if(isec==0) {
sensorCentroidX += xCheck;
sensorCentroidZ += zCheck;
sensorCount++;
};
// placement (note: transformations are in reverse order)
// - transformations operate on global coordinates; the corresponding
// generator coordinates are provided in the comments
auto sensorPV = gasvolVol.placeVolume(sensorVol,
RotationZ(sectorRotation) // rotate about beam axis to sector
* Translation3D(sensorSphPos.x(), sensorSphPos.y(), sensorSphPos.z()) // move sphere to reference position
* RotationX(phiGen) // rotate about `zGen`
* RotationZ(thetaGen) // rotate about `yGen`
* Translation3D(sensorSphRadius, 0., 0.) // push radially to spherical surface
* RotationY(M_PI/2) // rotate sensor to be compatible with generator coords
* RotationZ(-M_PI/2) // correction for readout segmentation mapping
);
// generate LUT for module number -> sensor position, for readout mapping tests
//if(isec==0) printf("%d %f %f\n",imod,sensorPV.position().x(),sensorPV.position().y());
// properties
sensorPV.addPhysVolID("sector", isec).addPhysVolID("module", imod);
DetElement sensorDE(det, Form("sensor_de%d_%d", isec, imod), (imod<<3)|isec ); // id must match IRTAlgorithm usage
sensorDE.setPlacement(sensorPV);
if(!debug_optics) {
SkinSurface sensorSkin(desc, sensorDE, Form("sensor_optical_surface%d", isec), sensorSurf, sensorVol);
sensorSkin.isValid();
};
// increment sensor module number
imod++;
}; // end patch cuts
}; // end phiGen loop
}; // end thetaGen loop
// calculate centroid sensor position
if(isec==0) {
sensorCentroidX /= sensorCount;
sensorCentroidZ /= sensorCount;
};
// END SENSOR MODULE LOOP ------------------------
// BUILD MIRRORS ====================================================================
// derive spherical mirror parameters `(zM,xM,rM)`, for given image point
// coordinates `(zI,xI)` and `dO`, defined as the z-distance between the
// object and the mirror surface
// - all coordinates are specified w.r.t. the object point coordinates
// - this is point-to-point focusing, but it can be used to effectively steer
// parallel-to-point focusing
double zM,xM,rM;
auto FocusMirror = [&zM,&xM,&rM](double zI, double xI, double dO) {
zM = dO*zI / (2*dO-zI);
xM = dO*xI / (2*dO-zI);
rM = dO - zM;
};
// attributes, re-defined w.r.t. IP, needed for mirror positioning
double zS = sensorSphCenterZ + vesselZmin; // sensor sphere attributes
double xS = sensorSphCenterX;
double rS = sensorSphRadius;
double B = vesselZmax - mirrorBackplane; // distance between IP and mirror back plane
// focus 1: set mirror to focus IP on center of sensor sphere `(zS,xS)`
/*double zF = zS;
double xF = xS;
FocusMirror(zF,xF,B);*/
// focus 2: move focal region along sensor sphere radius, according to `focusTuneLong`
// - specifically, along the radial line which passes through the approximate centroid
// of the sensor region `(sensorCentroidZ,sensorCentroidX)`
// - `focusTuneLong` is the distance to move, given as a fraction of `sensorSphRadius`
// - `focusTuneLong==0` means `(zF,xF)==(zS,xS)`
// - `focusTuneLong==1` means `(zF,xF)` will be on the sensor sphere, near the centroid
/*
double zC = sensorCentroidZ + vesselZmin;
double xC = sensorCentroidX;
double slopeF = (xC-xS) / (zC-zS);
double thetaF = std::atan(std::fabs(slopeF));
double zF = zS + focusTuneLong * sensorSphRadius * std::cos(thetaF);
double xF = xS - focusTuneLong * sensorSphRadius * std::sin(thetaF);
//FocusMirror(zF,xF,B);
// focus 3: move along line perpendicular to focus 2's radial line,
// according to `focusTunePerp`, with the same numerical scale as `focusTuneLong`
zF += focusTunePerp * sensorSphRadius * std::cos(M_PI/2-thetaF);
xF += focusTunePerp * sensorSphRadius * std::sin(M_PI/2-thetaF);
FocusMirror(zF,xF,B);
*/
// focus 4: use (z,x) coordinates for tune parameters
double zF = zS + focusTuneZ;
double xF = xS + focusTuneX;
FocusMirror(zF,xF,B);
// re-define mirror attributes to be w.r.t vessel front plane
double mirrorCenterZ = zM - vesselZmin;
double mirrorCenterX = xM;
double mirrorRadius = rM;
// spherical mirror patch cuts and rotation
double mirrorThetaRot = std::asin(mirrorCenterX/mirrorRadius);
double mirrorTheta1 = mirrorThetaRot - std::asin((mirrorCenterX-mirrorRmin)/mirrorRadius);
double mirrorTheta2 = mirrorThetaRot + std::asin((mirrorRmax-mirrorCenterX)/mirrorRadius);
// if debugging, draw full sphere
if(debug_mirror) { mirrorTheta1=0; mirrorTheta2=M_PI; /*mirrorPhiw=2*M_PI;*/ };
// solid : create sphere at origin, with specified angular limits;
// phi limits are increased to fill gaps (overlaps are cut away later)
Sphere mirrorSolid1(
mirrorRadius,
mirrorRadius + mirrorThickness,
mirrorTheta1,
mirrorTheta2,
-40*degree,
40*degree
);
/* CAUTION: if any of the relative placements or boolean operations below
* are changed, you MUST make sure this does not break access to the sphere
* primitive and positioning in Juggler `IRTAlgorithm`; cross check the
* mirror sphere attributes carefully!
*/
/*
// PRINT MIRROR ATTRIBUTES (before any sector z-rotation)
printf("zM = %f\n",zM); // sphere centerZ, w.r.t. IP
printf("xM = %f\n",xM); // sphere centerX, w.r.t. IP
printf("rM = %f\n",rM); // sphere radius
*/
// mirror placement transformation (note: transformations are in reverse order)
auto mirrorPos = Position(mirrorCenterX, 0., mirrorCenterZ) + originFront;
Transform3D mirrorPlacement(
Translation3D(mirrorPos.x(), mirrorPos.y(), mirrorPos.z()) // re-center to specified position
* RotationY(-mirrorThetaRot) // rotate about vertical axis, to be within vessel radial walls
);
// cut overlaps with other sectors using "pie slice" wedges, to the extent specified
// by `mirrorPhiw`
Tube pieSlice( 0.01*cm, vesselRmax2, tankLength/2.0, -mirrorPhiw/2.0, mirrorPhiw/2.0);
IntersectionSolid mirrorSolid2( pieSlice, mirrorSolid1, mirrorPlacement );
// mirror volume, attributes, and placement
Volume mirrorVol(detName+"_mirror_"+secName, mirrorSolid2, mirrorMat);
mirrorVol.setVisAttributes(mirrorVis);
auto mirrorPV2 = gasvolVol.placeVolume(mirrorVol,
RotationZ(sectorRotation) // rotate about beam axis to sector
* Translation3D(0,0,0)
);
// properties
DetElement mirrorDE(det, Form("mirror_de%d", isec), isec);
mirrorDE.setPlacement(mirrorPV2);
SkinSurface mirrorSkin(desc, mirrorDE, Form("mirror_optical_surface%d", isec), mirrorSurf, mirrorVol);
mirrorSkin.isValid();
}; // END SECTOR LOOP //////////////////////////
// place gas volume
PlacedVolume gasvolPV = vesselVol.placeVolume(gasvolVol,Position(0, 0, 0));
DetElement gasvolDE(det, "gasvol_de", 0);
gasvolDE.setPlacement(gasvolPV);
// place mother volume (vessel)
Volume motherVol = desc.pickMotherVolume(det);
PlacedVolume vesselPV = motherVol.placeVolume(vesselVol,
Position(0, 0, vesselZmin) - originFront
);
vesselPV.addPhysVolID("system", detID);
det.setPlacement(vesselPV);
return det;
};
// clang-format off
DECLARE_DETELEMENT(athena_DRICH, createDetector)
src/FieldMapBrBz.cpp 0 → 100644
View file @ 1008a184
#include <DD4hep/DetFactoryHelper.h>
#include <DD4hep/FieldTypes.h>
#include <DD4hep/Printout.h>
#include <XML/Utilities.h>
#include <cstdlib>
#include <filesystem>
#include <fstream>
#include <iostream>
#include <sstream>
#include <stdexcept>
#include <string>
#include <tuple>
namespace fs = std::filesystem;
#include "FileLoaderHelper.h"
using namespace dd4hep;
// implementation of the field map
class FieldMapBrBz : public dd4hep::CartesianField::Object
{
public:
FieldMapBrBz(const std::string &field_type = "magnetic");
void Configure(double rmin, double rmax, double rstep, double zmin, double zmax, double zstep);
void LoadMap(const std::string &map_file, double scale);
void GetIndices(double r, double z, int &ir, int &iz, double &dr, double &dz);
void SetTransform(const Transform3D &tr) { trans = tr; trans_inv = tr.Inverse(); }
virtual void fieldComponents(const double *pos, double *field);
private:
Transform3D trans, trans_inv;
double rmin, rmax, rstep, zmin, zmax, zstep;
std::vector<std::vector<std::vector<double>>> Bvals;
};
// constructor
FieldMapBrBz::FieldMapBrBz(const std::string &field_type)
{
std::string ftype = field_type;
for (auto &c : ftype) { c = tolower(c); }
// set type
if (ftype == "magnetic") {
type = CartesianField::MAGNETIC;
} else if (ftype == "electric") {
type = CartesianField::ELECTRIC;
} else {
type = CartesianField::UNKNOWN;
std::cout << "FieldMapBrBz Warning: Unknown field type " << field_type << "!" << std::endl;
}
}
void FieldMapBrBz::Configure(double r1, double r2, double rs, double z1, double z2, double zs)
{
rmin = r1;
rmax = r2;
rstep = rs;
zmin = z1;
zmax = z2;
zstep = zs;
int nr = int((r2 - r1)/rs) + 2;
int nz = int((z2 - z1)/zs) + 2;
Bvals.resize(nr);
for (auto &B2 : Bvals) {
B2.resize(nz);
for (auto &B : B2) {
B.resize(2, 0.);
}
}
}
void FieldMapBrBz::GetIndices(double r, double z, int &ir, int &iz, double &dr, double &dz)
{
// boundary check
if (r > rmax || r < rmin || z > zmax || z < zmin) {
ir = -1;
iz = -1;
return;
}
// get indices
double idr, idz;
dr = std::modf((r - rmin)/rstep, &idr);
dz = std::modf((z - zmin)/zstep, &idz);
ir = static_cast<int>(idr);
iz = static_cast<int>(idz);
}
// load data
void FieldMapBrBz::LoadMap(const std::string &map_file, double scale)
{
std::string line;
std::ifstream input(map_file);
if (!input) {
std::cout << "FieldMapBrBz Error: file \"" << map_file << "\" cannot be read." << std::endl;
}
double r, z, br, bz;
int ir, iz;
double dr, dz;
while (std::getline(input, line).good()) {
std::istringstream iss(line);
iss >> r >> z >> br >> bz;
GetIndices(r, z, ir, iz, dr, dz);
if (ir < 0 || iz < 0) {
std::cout << "FieldMapBrBz Warning: coordinates out of range ("
<< r << ", " << z << "), skipped it." << std::endl;
} else {
Bvals[ir][iz] = {br*scale, bz*scale};
// ROOT::Math::XYZPoint p(r, 0, z);
// std::cout << p << " -> " << trans*p << std::endl;
// std::cout << ir << ", " << iz << ", " << br << ", " << bz << std::endl;
}
}
}
// get field components
void FieldMapBrBz::fieldComponents(const double *pos, double *field)
{
// coordinate conversion
auto p = trans_inv*ROOT::Math::XYZPoint(pos[0], pos[1], pos[2]);
// coordinates conversion
const double r = sqrt(p.x()*p.x() + p.y()*p.y());
const double z = p.z();
const double phi = atan2(p.y(), p.x());
int ir, iz;
double dr, dz;
GetIndices(r, z, ir, iz, dr, dz);
// out of the range
if (ir < 0 || iz < 0) { return; }
// p1 p3
// p
// p0 p2
auto &p0 = Bvals[ir][iz];
auto &p1 = Bvals[ir][iz + 1];
auto &p2 = Bvals[ir + 1][iz];
auto &p3 = Bvals[ir + 1][iz + 1];
// linear interpolation
double Br = p0[0] * (1-dr) * (1-dz)
+ p1[0] * (1-dr) * dz
+ p2[0] * dr * (1-dz)
+ p3[0] * dr * dz;
double Bz = p0[1] * (1-dr) * (1-dz)
+ p1[1] * (1-dr) * dz
+ p2[1] * dr * (1-dz)
+ p3[1] * dr * dz;
// convert Br Bz to Bx By Bz
auto B = trans*ROOT::Math::XYZPoint(Br*sin(phi), Br*cos(phi), Bz);
field[0] += B.x()*tesla;
field[1] += B.y()*tesla;
field[2] += B.z()*tesla;
return;
}
// assign the field map to CartesianField
static Ref_t create_field_map_brbz(Detector & /*lcdd*/, xml::Handle_t handle)
{
xml_comp_t x_par(handle);
if (!x_par.hasAttr(_Unicode(field_map))) {
throw std::runtime_error("FieldMapBrBz Error: must have an xml attribute \"field_map\" for the field map.");
}
CartesianField field;
std::string field_type = x_par.attr<std::string>(_Unicode(field_type));
// dimensions
xml_comp_t x_dim = x_par.dimensions();
// min, max, step
xml_comp_t r_dim = x_dim.child(_Unicode(transverse));
xml_comp_t z_dim = x_dim.child(_Unicode(longitudinal));
std::string field_map_file = x_par.attr<std::string>(_Unicode(field_map));
std::string field_map_url = x_par.attr<std::string>(_Unicode(url));
EnsureFileFromURLExists(field_map_url,field_map_file);
double field_map_scale = x_par.attr<double>(_Unicode(scale));
if( !fs::exists(fs::path(field_map_file)) ) {
printout(ERROR, "FieldMapBrBz", "file " + field_map_file + " does not exist");
printout(ERROR, "FieldMapBrBz", "use a FileLoader plugin before the field element");
std::_Exit(EXIT_FAILURE);
}
auto map = new FieldMapBrBz(field_type);
map->Configure(r_dim.rmin(), r_dim.rmax(), r_dim.step(), z_dim.zmin(), z_dim.zmax(), z_dim.step());
// translation, rotation
static double deg2r = ROOT::Math::Pi()/180.;
RotationZYX rot(0., 0., 0.);
if (x_dim.hasChild(_Unicode(rotation))) {
xml_comp_t rot_dim = x_dim.child(_Unicode(rotation));
rot = RotationZYX(rot_dim.z()*deg2r, rot_dim.y()*deg2r, rot_dim.x()*deg2r);
}
Translation3D trans(0., 0., 0.);
if (x_dim.hasChild(_Unicode(translation))) {
xml_comp_t trans_dim = x_dim.child(_Unicode(translation));
trans = Translation3D(trans_dim.x(), trans_dim.y(), trans_dim.z());
}
map->SetTransform(trans*rot);
map->LoadMap(field_map_file, field_map_scale);
field.assign(map, x_par.nameStr(), "FieldMapBrBz");
return field;
}
DECLARE_XMLELEMENT(FieldMapBrBz, create_field_map_brbz)
src/FileLoader.cpp 0 → 100644
View file @ 1008a184
#include <DD4hep/DetFactoryHelper.h>
#include <DD4hep/Primitives.h>
#include <DD4hep/Factories.h>
#include <DD4hep/Printout.h>
#include <XML/Utilities.h>
#include <fmt/core.h>
#include <filesystem>
#include <iostream>
#include <cstdlib>
#include <string>
#include "FileLoaderHelper.h"
using namespace dd4hep;
void usage(int argc, char** argv) {
std::cout <<
"Usage: -plugin <name> -arg [-arg] \n"
" name: factory name FileLoader \n"
" cache:<string> cache location (may be read-only) \n"
" file:<string> file location \n"
" url:<string> url location \n"
" cmd:<string> download command with {0} for url, {1} for output \n"
"\tArguments given: " << arguments(argc,argv) << std::endl;
std::exit(EINVAL);
}
// Plugin to download files
long load_file(
Detector& /* desc */,
int argc,
char** argv
) {
// argument parsing
std::string cache, file, url;
std::string cmd("curl --retry 5 -f {0} -o {1}");
for (int i = 0; i < argc && argv[i]; ++i) {
if (0 == std::strncmp("cache:", argv[i], 6)) cache = (argv[i] + 6);
else if (0 == std::strncmp("file:", argv[i], 5)) file = (argv[i] + 5);
else if (0 == std::strncmp("url:", argv[i], 4)) url = (argv[i] + 4);
else if (0 == std::strncmp("cmd:", argv[i], 4)) cmd = (argv[i] + 4);
else usage(argc, argv);
}
printout(DEBUG, "FileLoader", "arg cache: " + cache);
printout(DEBUG, "FileLoader", "arg file: " + file);
printout(DEBUG, "FileLoader", "arg url: " + url);
printout(DEBUG, "FileLoader", "arg cmd: " + cmd);
// if file or url is empty, do nothing
if (file.empty()) {
printout(WARNING, "FileLoader", "no file specified");
return 0;
}
if (url.empty()) {
printout(WARNING, "FileLoader", "no url specified");
return 0;
}
EnsureFileFromURLExists(url, file, cache, cmd);
return 0;
}
DECLARE_APPLY(FileLoader, load_file)
src/FileLoaderHelper.h 0 → 100644
View file @ 1008a184
#pragma once
#include <DD4hep/DetFactoryHelper.h>
#include <DD4hep/Primitives.h>
#include <DD4hep/Factories.h>
#include <DD4hep/Printout.h>
#include <fmt/core.h>
#include <filesystem>
#include <iostream>
#include <cstdlib>
#include <string>
namespace fs = std::filesystem;
using namespace dd4hep;
// Function to download files
inline void
EnsureFileFromURLExists(
std::string url,
std::string file,
std::string cache = "",
std::string cmd = "curl --retry 5 -f {0} -o {1}"
) {
// parse cache for environment variables
auto pos = std::string::npos;
while ((pos = cache.find('$')) != std::string::npos) {
auto after = cache.find_first_not_of(
"ABCDEFGHIJKLMNOPQRSTUVWXYZ"
"abcdefghijklmnopqrstuvwxyz"
"0123456789"
"_",
pos + 1);
if (after == std::string::npos) after = cache.size(); // cache ends on env var
const std::string env_name(cache.substr(pos + 1, after - pos - 1));
auto env_ptr = std::getenv(env_name.c_str());
const std::string env_value(env_ptr != nullptr ? env_ptr : "");
cache.erase(pos, after - pos);
cache.insert(pos, env_value);
printout(INFO, "FileLoader", "$" + env_name + " -> " + env_value);
}
// create file path
fs::path file_path(file);
// create hash from url, hex of unsigned long long
std::string hash = fmt::format("{:016x}", dd4hep::detail::hash64(url)); // TODO: Use c++20 std::fmt
// create file parent path, if not exists
fs::path parent_path = file_path.parent_path();
if (!fs::exists(parent_path)) {
if (fs::create_directories(parent_path) == false) {
printout(ERROR, "FileLoader", "parent path " + parent_path.string() + " cannot be created");
printout(ERROR, "FileLoader", "check permissions and retry");
std::_Exit(EXIT_FAILURE);
}
}
// if file exists and is symlink to correct hash
fs::path hash_path(parent_path / hash);
if (fs::exists(file_path)
&& fs::equivalent(file_path, hash_path)) {
printout(INFO, "FileLoader", "Link " + file + " -> hash " + hash + " already exists");
return;
}
// if hash does not exist, we try to retrieve file from cache
if (!fs::exists(hash_path)) {
// recursive loop into cache directory
fs::path cache_path(cache);
printout(INFO, "FileLoader", "Cache " + cache_path.string());
if (fs::exists(cache_path)) {
for (auto const& dir_entry: fs::recursive_directory_iterator(cache_path)) {
if (!dir_entry.is_directory()) continue;
fs::path cache_dir_path = cache_path / dir_entry;
printout(INFO, "FileLoader", "Checking " + cache_dir_path.string());
fs::path cache_hash_path = cache_dir_path / hash;
if (fs::exists(cache_hash_path)) {
// symlink hash to cache/.../hash
printout(INFO, "FileLoader", "File " + file + " with hash " + hash + " found in " + cache_hash_path.string());
try {
fs::create_symlink(cache_hash_path, hash_path);
} catch (const fs::filesystem_error&) {
printout(ERROR, "FileLoader", "unable to link from " + hash_path.string() + " to " + cache_hash_path.string());
printout(ERROR, "FileLoader", "check permissions and retry");
std::_Exit(EXIT_FAILURE);
}
break;
}
}
}
}
// if hash does not exist, we try to retrieve file from url
if (!fs::exists(hash_path)) {
cmd = fmt::format(cmd, url, hash_path.c_str()); // TODO: Use c++20 std::fmt
printout(INFO, "FileLoader", "Downloading " + file + " as hash " + hash + " with " + cmd);
// run cmd
auto ret = std::system(cmd.c_str());
if (!fs::exists(hash_path)) {
printout(ERROR, "FileLoader", "unable to run cmd " + cmd);
printout(ERROR, "FileLoader", "check command and retry");
printout(ERROR, "FileLoader", "hint: allow insecure connections with -k");
std::_Exit(EXIT_FAILURE);
}
}
// check if file already exists
if (fs::exists(file_path)) {
// file already exists
if (fs::is_symlink(file_path)) {
// file is symlink
if (fs::equivalent(hash_path, fs::read_symlink(file_path))) {
// link points to correct path
return;
} else {
// link points to incorrect path
if (fs::remove(file_path) == false) {
printout(ERROR, "FileLoader", "unable to remove symlink " + file_path.string());
printout(ERROR, "FileLoader", "check permissions or remove manually");
std::_Exit(EXIT_FAILURE);
}
}
} else {
// file exists but not symlink
printout(ERROR, "FileLoader", "will not remove actual file " + file_path.string());
printout(ERROR, "FileLoader", "check content, remove manually, and retry");
std::_Exit(EXIT_FAILURE);
}
}
// file_path now does not exist
// symlink file_path to hash_path
try {
// use new path from hash so file link is local
fs::create_symlink(fs::path(hash), file_path);
} catch (const fs::filesystem_error&) {
printout(ERROR, "FileLoader", "unable to link from " + file_path.string() + " to " + hash_path.string());
printout(ERROR, "FileLoader", "check permissions and retry");
std::_Exit(EXIT_FAILURE);
}
}
src/ForwardRICH_geo.cpp deleted 100644 → 0
View file @ f87c4897
//==========================================================================
// Forward Ring Imaging Cherenkov Detector
//--------------------------------------------------------------------------
//
// Author: C. Peng (ANL)
// Date: 09/30/2020
//
//==========================================================================
#include <XML/Helper.h>
#include "TMath.h"
#include "TString.h"
#include "ref_utils.h"
#include "Math/Point2D.h"
#include "DDRec/Surface.h"
#include "DDRec/DetectorData.h"
#include "DD4hep/OpticalSurfaces.h"
#include "DD4hep/DetFactoryHelper.h"
#include "DD4hep/Printout.h"
using namespace std;
using namespace dd4hep;
using namespace dd4hep::rec;
// create the detector
static Ref_t createDetector(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);
xml::Component dims = detElem.dimensions();
xml::Component rads = detElem.child(_Unicode(radiator));
xml::Component mir = detElem.child(_Unicode(mirror));
xml::Component mcp = detElem.child(_Unicode(mcppmt));
xml::Component tank = detElem.child(_Unicode(tank));
// dimensions
double z0 = dims.z0();
double length = dims.length();
double rmin = dims.rmin();
double rmax0 = dims.attr<double>(_Unicode(rmax0));
double rmax1 = dims.attr<double>(_Unicode(rmax1));
double rmax2 = dims.attr<double>(_Unicode(rmax2));
double snout_length = dims.attr<double>(_Unicode(snout_length));
// mirror setting
auto mThick = mir.thickness();
auto mirZ = mir.attr<double>(_Unicode(z0));
// mcppmt setting
auto pRmin = mcp.rmin();
auto pRmax = mcp.rmax();
auto pThick = mcp.thickness();
auto pSize = mcp.attr<double>(_Unicode(module_size));
auto pGap = mcp.attr<double>(_Unicode(module_gap));
auto pTol = mcp.attr<double>(_Unicode(rtol));
auto pZ = mcp.attr<double>(_Unicode(z0));
// tank parameters
double tank_length = length - snout_length;
// materials
auto mirMat = desc.material(mir.materialStr());
auto gasMat = desc.material(rads.materialStr());
auto mcpMat = desc.material(mcp.materialStr());
double front_offset = snout_length+tank_length/2.0;
// constants
auto richCenterAngle = std::atan((rmin + (rmax1 - rmin)/2.)/(front_offset+mirZ));
//std::cout << richCenterAngle*180./M_PI << std::endl;
// an envelope for the detector
// use a complicated shape to avoid conflict with the other parts
// cone for radiator and the first set of mirrors
double halfLength = length/2.;
Cone env1(snout_length/2.0, rmin, rmax0, rmin, rmax1);
// envelope for detection plane
// Cone env2(halfLength - pZ/2., rmin, pRmax, rmin, rmax2);
Tube env2(rmin, pRmax + pTol + pGap + 1.0*cm, tank_length/2., 0., 2*M_PI);
UnionSolid envShape(env2, env1, Position(0., 0., -tank_length/2.-snout_length/2));
Volume envVol(detName + "_envelope", envShape, gasMat);
envVol.setVisAttributes(desc.visAttributes(detElem.visStr()));
// ---------------
// spherical mirrors inside it
int ilayer = 1;
// optical surface
OpticalSurfaceManager surfMgr = desc.surfaceManager();
OpticalSurface mirSurf = surfMgr.opticalSurface("MirrorOpticalSurface");
// mirror slices
int imod = 1;
for (xml::Collection_t sl(mir, _Unicode(slice)); sl; ++sl, ++imod) {
auto focus = sl.attr<double>(_Unicode(focus));
auto wphi = sl.attr<double>(_Unicode(phiw));
auto rotZ = sl.attr<double>(_Unicode(rotz));
auto mRmin = sl.attr<double>(_Unicode(rmin));
auto mRmax = sl.attr<double>(_Unicode(rmax));
double curve = 0.;
if (sl.hasAttr(_Unicode(curve))) {
curve = sl.attr<double>(_Unicode(curve));
}
// geometry of mirror slice
PlacedVolume mirPV;
Volume mirVol(Form("mirror_v_dummy%d", imod));
mirVol.setMaterial(mirMat);
mirVol.setVisAttributes(desc.visAttributes(mir.visStr()));
// spherical mirror
if (curve > 0.) {
// somehow geant4 does not support -wphi/2. to wphi/2., so additonal rotation in Z
double mTheta1 = std::asin(mRmin/curve);
double mTheta2 = std::asin(mRmax/curve);
double rotY = -std::asin(focus/curve);
mirVol.setSolid(Sphere(curve, curve + mThick, mTheta1, mTheta2, 0., wphi));
// action is in a reverse order
Transform3D tr = Translation3D(0., 0., mirZ - front_offset) // move for z position
* RotationZ(rotZ) // rotate phi angle
* RotationY(rotY) // rotate for focus point
* RotationX(180*degree)
* Translation3D(0., 0., -curve) // move spherical shell to origin
* RotationZ(-wphi/2.); // center phi angle to 0. (-wphi/2., wphi/2.)
mirPV = envVol.placeVolume(mirVol, tr);
// plane mirror
} else {
mirVol.setSolid(Tube(mRmin, mRmax, mThick/2.0, 0., wphi));
Transform3D tr = Translation3D(0., 0., mirZ - front_offset) // move for z position
* RotationZ(rotZ) // rotate phi angle
* RotationZ(-wphi/2.); // center phi angle to 0. (-wphi/2., wphi/2.)
mirPV = envVol.placeVolume(mirVol, tr);
}
mirPV.addPhysVolID("layer", ilayer).addPhysVolID("module", imod);
DetElement mirDE(det, Form("Mirror_DE%d", imod), imod);
mirDE.setPlacement(mirPV);
SkinSurface mirSurfBorder(desc, mirDE, Form("RICHmirror%d", imod), mirSurf, mirVol);
mirSurfBorder.isValid();
}
ilayer++;
// ---------------
// photo-detector unit
// Fill the photo-detection plane with square shape MCP-PMTs
Box mcpShape1(pSize/2.0, pSize/2.0, pThick/2.0);
Volume mcpVol1("mcppmt_v_material", mcpShape1, mcpMat);
// a thin layer of cherenkov gas for accepting optical photons
Box mcpShape(pSize/2.0, pSize/2.0, pThick/2.0 + 0.1*mm);
Volume mcpVol("mcppmt_v", mcpShape, gasMat);
mcpVol.placeVolume(mcpVol1, Position(0., 0., -0.1*mm));
mcpVol.setVisAttributes(desc.visAttributes(mcp.visStr()));
sens.setType("photoncounter");
mcpVol.setSensitiveDetector(sens);
// photo-detector plane envelope
for (size_t ipd = 0; ipd < 6; ++ipd) {
double phmin = -M_PI/6.5; // added 0.5 to make it smaller
double phmax = M_PI/6.5;
Tube pdEnvShape(pRmin - pTol - pGap, pRmax + pTol + pGap, pThick/2.0 + 0.1*cm, phmin, phmax);
Volume pdVol("pd_envelope", pdEnvShape, desc.material("AirOptical"));
auto points = ref::utils::fillSquares({0., 0.}, pSize + pGap, pRmin + pTol + pGap, pRmax + pTol + pGap, phmin, phmax);
for (size_t i = 0; i < points.size(); ++i) {
auto pt = points[i];
auto mcpPV = pdVol.placeVolume(mcpVol, Position(pt.x(), pt.y(), 0.));
mcpPV.addPhysVolID("layer", ilayer).addPhysVolID("module", i + 1);
DetElement mcpDE(det, Form("MCPPMT_DE%d_%d", ipd + 1, i + 1), i + 1);
mcpDE.setPlacement(mcpPV);
}
Transform3D tr = Translation3D(0., 0., -front_offset + pZ + pThick/2.0) // move for z position
* RotationZ(ipd*M_PI/3.) // rotate phi angle
* RotationY(-richCenterAngle); // rotate to perpendicular position
auto pdPV = envVol.placeVolume(pdVol, tr);
pdPV.addPhysVolID("layer", ilayer).addPhysVolID("piece", ipd + 1);
}
Volume motherVol = desc.pickMotherVolume(det);
PlacedVolume envPV = motherVol.placeVolume(envVol, Position(0, 0, z0 + halfLength));
envPV.addPhysVolID("system", detID);
det.setPlacement(envPV);
return det;
}
//@}
// clang-format off
DECLARE_DETELEMENT(refdet_ForwardRICH, createDetector)
src/GaseousRICH_geo.cpp 0 → 100644
View file @ 1008a184
//==========================================================================
// Gaseous Ring Imaging Cherenkov Detector
//--------------------------------------------------------------------------
//
// Author: C. Peng (ANL)
// Date: 09/30/2020
//
//==========================================================================
#include <XML/Helper.h>
#include "TMath.h"
#include "TString.h"
#include "GeometryHelpers.h"
#include "Math/Point2D.h"
#include "DDRec/Surface.h"
#include "DDRec/DetectorData.h"
#include "DD4hep/OpticalSurfaces.h"
#include "DD4hep/DetFactoryHelper.h"
#include "DD4hep/Printout.h"
using namespace std;
using namespace dd4hep;
using namespace dd4hep::rec;
// headers
void build_radiator(Detector &desc, Volume &env, xml::Component plm, const Position &offset);
void build_mirrors(Detector &desc, DetElement &sdet, Volume &env, xml::Component plm, const Position &offset);
void build_sensors(Detector &desc, Volume &env, xml::Component plm, const Position &offset, SensitiveDetector &sens);
// 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;
}
// create the detector
static Ref_t createDetector(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);
xml::Component dims = detElem.dimensions();
// build a big envelope for all the components, filled with optical material to ensure transport of optical photons
double z0 = dims.z0();
double length = dims.length();
double rmin = dims.rmin();
double rmax0 = dims.attr<double>(_Unicode(rmax0));
double rmax1 = dims.attr<double>(_Unicode(rmax1));
double rmax2 = dims.attr<double>(_Unicode(rmax2));
double snout_length = dims.attr<double>(_Unicode(snout_length));
// fill envelope with radiator materials (need optical property for optical photons)
auto gasMat = desc.material(dd4hep::getAttrOrDefault(detElem, _Unicode(gas), "AirOptical"));
Cone snout(snout_length/2.0, rmin, rmax0, rmin, rmax1);
Tube tank(rmin, rmax2, (length - snout_length)/2., 0., 2*M_PI);
// shift the snout to the left side of tank
UnionSolid envShape(tank, snout, Position(0., 0., -length/2.));
// some reference points
auto snout_front = Position(0., 0., -(length + snout_length)/2.);
// auto snout_end = Position(0., 0., -(length - snout_length)/2.); // tank_front
// auto tank_end = Position(0., 0., (length - snout_length)/2.);
Volume envVol(detName + "_envelope", envShape, gasMat);
envVol.setVisAttributes(desc.visAttributes(detElem.visStr()));
// sensitive detector type
sens.setType("tracker");
// @TODO: place a radiator
// build_radiator(desc, envVol, detElement.child(_Unicode(radiator)), snout_front);
// place mirrors
build_mirrors(desc, det, envVol, detElem.child(_Unicode(mirrors)), snout_front);
// place photo-sensors
build_sensors(desc, envVol, detElem.child(_Unicode(sensors)), snout_front, sens);
Volume motherVol = desc.pickMotherVolume(det);
PlacedVolume envPV = motherVol.placeVolume(envVol, Position(0, 0, z0) - snout_front);
envPV.addPhysVolID("system", detID);
det.setPlacement(envPV);
return det;
}
// @TODO: implement a radiator, now the envelope serves as the radiator
void build_radiator(Detector &desc, Volume &env, xml::Component plm, const Position &offset)
{
// place holder
}
// place mirrors
void build_mirrors(Detector &desc, DetElement &sdet, Volume &env, xml::Component plm, const Position &offset)
{
double thickness = dd4hep::getAttrOrDefault<double>(plm, _Unicode(dz), 1.*dd4hep::mm);
auto mat = desc.material(plm.attr<std::string>(_Unicode(material)));
auto vis = desc.visAttributes(plm.attr<std::string>(_Unicode(vis)));
// optical surface
OpticalSurfaceManager surfMgr = desc.surfaceManager();
auto surf = surfMgr.opticalSurface(dd4hep::getAttrOrDefault(plm, _Unicode(surface), "MirrorOpticalSurface"));
// placements
auto gpos = get_xml_xyz(plm, _Unicode(position)) + offset;
auto grot = get_xml_xyz(plm, _Unicode(position));
int imod = 1;
for (xml::Collection_t mir(plm, _Unicode(mirror)); mir; ++mir, ++imod) {
double rmin = mir.attr<double>(_Unicode(rmin));
double rmax = mir.attr<double>(_Unicode(rmax));
double phiw = dd4hep::getAttrOrDefault<double>(mir, _Unicode(phiw), 2.*M_PI);
double rotz = dd4hep::getAttrOrDefault<double>(mir, _Unicode(rotz), 0.);
double roty = dd4hep::getAttrOrDefault<double>(mir, _Unicode(roty), 0.);
double rotx = dd4hep::getAttrOrDefault<double>(mir, _Unicode(rotx), 0.);
Volume vol(Form("mirror_v_dummy%d", imod));
vol.setMaterial(mat);
vol.setVisAttributes(vis);
// mirror curvature
double curve = dd4hep::getAttrOrDefault<double>(mir, _Unicode(curve), 0.);
// spherical mirror
if (curve > 0.) {
double th1 = std::asin(rmin/curve);
double th2 = std::asin(rmax/curve);
vol.setSolid(Sphere(curve, curve + thickness, th1, th2, 0., phiw));
// plane mirror
} else {
vol.setSolid(Tube(rmin, rmax, thickness/2., 0., phiw));
}
// transforms are in a reverse order
Transform3D tr = Translation3D(gpos.x(), gpos.y(), gpos.z())
* RotationZYX(grot.z(), grot.y(), grot.x())
* RotationZ(rotz) // rotation of the piece
* RotationY(roty) // rotation of the piece
* RotationX(rotx) // rotation of the piece
* Translation3D(0., 0., -curve) // move spherical shell to origin (no move for planes)
* RotationZ(-phiw/2.); // center phi angle to 0. (-phiw/2., phiw/2.)
auto pv = env.placeVolume(vol, tr);
DetElement de(sdet, Form("mirror_de%d", imod), imod);
de.setPlacement(pv);
SkinSurface skin(desc, de, Form("mirror_optical_surface%d", imod), surf, vol);
skin.isValid();
}
}
// place photo-sensors
void build_sensors(Detector &desc, Volume &env, xml::Component plm, const Position &offset, SensitiveDetector &sens)
{
// build sensor unit geometry
auto mod = plm.child(_Unicode(module));
double sx = mod.attr<double>(_Unicode(sx));
double sy = mod.attr<double>(_Unicode(sy));
double sz = mod.attr<double>(_Unicode(sz));
double gap = mod.attr<double>(_Unicode(gap));
auto mat = desc.material(mod.attr<std::string>(_Unicode(material)));
auto vis = desc.visAttributes(mod.attr<std::string>(_Unicode(vis)));
Box sensor(sx/2., sy/2., sz/2.);
Volume svol("sensor_v", sensor, mat);
svol.setVisAttributes(vis);
// a thin layer of cherenkov gas for accepting optical photons
auto opt = plm.child(_Unicode(optical));
double opthick = opt.attr<double>(_Unicode(thickness));
auto opmat = desc.material(opt.attr<std::string>(_Unicode(material)));
Box opshape(sx/2., sy/2., sz/2. + opthick/2.);
Volume opvol("sensor_v_optical", opshape, opmat);
opvol.placeVolume(svol, Position(0., 0., 0.));
opvol.setSensitiveDetector(sens);
// photo-detector plane envelope
auto gpos = get_xml_xyz(plm, _Unicode(position)) + offset;
auto grot = get_xml_xyz(plm, _Unicode(position));
int isec = 1;
for (xml::Collection_t sec(plm, _Unicode(sector)); sec; ++sec, ++isec) {
double rmin = sec.attr<double>(_Unicode(rmin));
double rmax = sec.attr<double>(_Unicode(rmax));
double phiw = dd4hep::getAttrOrDefault<double>(sec, _Unicode(phiw), 2.*M_PI);
double rotz = dd4hep::getAttrOrDefault<double>(sec, _Unicode(rotz), 0.);
double roty = dd4hep::getAttrOrDefault<double>(sec, _Unicode(roty), 0.);
double rotx = dd4hep::getAttrOrDefault<double>(sec, _Unicode(rotx), 0.);
// fill sensors to the piece
auto points = athena::geo::fillRectangles({0., 0.}, sx + gap, sy + gap, rmin - gap, rmax + gap, -phiw/2., phiw/2.);
int imod = 1;
for (auto &p : points) {
// transofrms are in a reversed order
Transform3D tr = Translation3D(gpos.x(), gpos.y(), gpos.z())
* RotationZYX(grot.z(), grot.y(), grot.x())
* RotationZ(rotz) // rotation of the sector
* RotationY(roty) // rotation of the sector
* RotationX(rotx) // rotation of the sector
* Translation3D(p.x(), p.y(), 0.); // place modules in each sector
auto pv = env.placeVolume(opvol, tr);
pv.addPhysVolID("sector", isec).addPhysVolID("module", imod++);
}
}
}
//@}
// clang-format off
DECLARE_DETELEMENT(athena_GaseousRICH, createDetector)
src/GeometryHelpers.cpp 0 → 100644
View file @ 1008a184
#include "GeometryHelpers.h"
// some utility functions that can be shared
namespace athena::geo {
typedef ROOT::Math::XYPoint Point;
// check if a 2d point is already in the container
bool already_placed(const Point &p, const std::vector<Point> &vec, double xs = 1.0, double ys = 1.0, double tol = 1e-6)
{
for (auto &pt : vec) {
if ((std::abs(pt.x() - p.x())/xs < tol) && std::abs(pt.y() - p.y())/ys < tol) {
return true;
}
}
return false;
}
// check if a square in a ring
inline bool rec_in_ring(const Point &pt, double sx, double sy, double rmin, double rmax, double phmin, double phmax)
{
if (pt.r() > rmax || pt.r() < rmin) {
return false;
}
// check four corners
std::vector<Point> pts {
Point(pt.x() - sx/2., pt.y() - sy/2.),
Point(pt.x() - sx/2., pt.y() + sy/2.),
Point(pt.x() + sx/2., pt.y() - sy/2.),
Point(pt.x() + sx/2., pt.y() + sy/2.),
};
for (auto &p : pts) {
if (p.r() > rmax || p.r() < rmin || p.phi() > phmax || p.phi() < phmin) {
return false;
}
}
return true;
}
// a helper function to recursively fill square in a ring
void add_rectangle(Point p, std::vector<Point> &res, double sx, double sy,
double rmin, double rmax, double phmin, double phmax, int max_depth = 20, int depth = 0)
{
// std::cout << depth << "/" << max_depth << std::endl;
// exceeds the maximum depth in searching or already placed
if ((depth > max_depth) || (already_placed(p, res, sx, sy))) {
return;
}
bool in_ring = rec_in_ring(p, sx, sy, rmin, rmax, phmin, phmax);
if (in_ring) {
res.emplace_back(p);
}
// continue search for a good placement or if no placement found yet
if (in_ring || res.empty()) {
// check adjacent squares
add_rectangle(Point(p.x() + sx, p.y()), res, sx, sy, rmin, rmax, phmin, phmax, max_depth, depth + 1);
add_rectangle(Point(p.x() - sx, p.y()), res, sx, sy, rmin, rmax, phmin, phmax, max_depth, depth + 1);
add_rectangle(Point(p.x(), p.y() + sy), res, sx, sy, rmin, rmax, phmin, phmax, max_depth, depth + 1);
add_rectangle(Point(p.x(), p.y() - sy), res, sx, sy, rmin, rmax, phmin, phmax, max_depth, depth + 1);
}
}
// fill squares
std::vector<Point> fillRectangles(Point ref, double sx, double sy, double rmin, double rmax, double phmin, double phmax)
{
// convert (0, 2pi) to (-pi, pi)
if (phmax > M_PI) {
phmin -= M_PI;
phmax -= M_PI;
}
// start with a seed square and find one in the ring
// move to center
ref = ref - Point(int(ref.x()/sx)*sx, int(ref.y()/sy)*sy);
std::vector<Point> res;
add_rectangle(ref, res, sx, sy, rmin, rmax, phmin, phmax, (int(rmax/sx) + 1)*(int(rmax/sy) + 1)*2);
return res;
}
// check if a regular polygon is inside a ring
bool poly_in_ring(const Point &p, int nsides, double lside, double rmin, double rmax, double phmin, double phmax)
{
// outer radius is contained
if ((p.r() + lside <= rmax) && (p.r() - lside >= rmin)) {
return true;
}
// inner radius is not contained
double rin = std::cos(M_PI/nsides)*lside;
if ((p.r() + rin > rmax) || (p.r() - rin < rmin)) {
return false;
}
// in between, check every corner
for (int i = 0; i < nsides; ++i) {
double phi = (i + 0.5)*2.*M_PI/static_cast<double>(nsides);
Point p2(p.x() + 2.*lside*std::sin(phi), p.y() + 2.*lside*std::cos(phi));
if ((p2.r() > rmax) || (p2.r() < rmin)) {
return false;
}
}
return true;
}
// recursively fill square (nside=4) or hexagon (nside=6) in a ring, other polygons won't work
void add_poly(Point p, std::vector<Point> &res, int nsides, double lside,
double rmin, double rmax, double phmin, double phmax, int max_depth = 20, int depth = 0)
{
// std::cout << depth << "/" << max_depth << std::endl;
// exceeds the maximum depth in searching or already placed
if ((depth > max_depth) || (already_placed(p, res, lside, lside))) {
return;
}
bool in_ring = poly_in_ring(p, nsides, lside, rmin, rmax, phmin, phmax);
if (in_ring) {
res.emplace_back(p);
}
// recursively add neigbors, continue if it was a good placement or no placement found yet
if (in_ring || res.empty()) {
double d = 2.*std::cos(M_PI/static_cast<double>(nsides))*lside;
for (int i = 0; i < nsides; ++i) {
double phi = i*2.*M_PI/static_cast<double>(nsides);
add_poly(Point(p.x() + 2.*lside*std::sin(phi), p.y() + 2.*lside*std::cos(phi)), res, nsides, lside,
rmin, rmax, phmin, phmax, max_depth, depth + 1);
}
}
}
std::vector<Point> fillHexagons(Point ref, double lside, double rmin, double rmax, double phmin, double phmax)
{
// convert (0, 2pi) to (-pi, pi)
if (phmax > M_PI) {
phmin -= M_PI;
phmax -= M_PI;
}
// start with a seed and find one in the ring
// move to center
ref = ref - Point(int(ref.x()/lside)*lside, int(ref.y()/lside)*lside);
std::vector<Point> res;
add_poly(ref, res, 6, lside, rmin, rmax, phmin, phmax, std::pow(int(rmax/lside) + 1, 2)*2);
return res;
}
} // athena::geo
src/GeometryHelpers.h 0 → 100644
View file @ 1008a184
#pragma once
#include <vector>
#include "Math/Point2D.h"
// some utility functions that can be shared
namespace athena::geo {
using Point = ROOT::Math::XYPoint;
/** Fill rectangles in a ring (disk).
*
* @param ref 2D reference point.
* @param sx x side length
* @param sy y side length
* @param rmin inner radius of disk
* @param rmax outer radius of disk to fill
* @param phmin phi min
* @param phmax phi max
*/
std::vector<Point> fillRectangles(Point ref, double sx, double sy, double rmin, double rmax,
double phmin = -M_PI, double phmax = M_PI);
// fill squares in a ring
inline std::vector<Point> fillSquares(Point ref, double size, double rmin, double rmax,
double phmin = -M_PI, double phmax = M_PI)
{
return fillRectangles(ref, size, size, rmin, rmax, phmin, phmax);
}
std::vector<Point> fillHexagons(Point ref, double lside, double rmin, double rmax,
double phmin = -M_PI, double phmax = M_PI);
} // athena::geo
src/HomogeneousCalorimeter_geo.cpp 0 → 100644
View file @ 1008a184
//==========================================================================
// 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>
using namespace dd4hep;
/** \addtogroup calorimeters Calorimeters
*/
/** \addtogroup Homogeneous Calorimeter
* \brief Type: **HomogeneousCalorimeter**.
* \author C. Peng
* \ingroup calorimeters
*
*
* \code
* <detector id="1" name="HyCal" type="HomogeneousCalorimeter" readout="EcalHits">
* <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" axis="x" begin="8" nmods="16"/>
* <line x="10.25*mm" y="30.75*mm" axis="y" begin="8" nmods="16"/>
* <line x="10.25*mm" y="51.25*mm" axis="z" begin="8" nmods="16"/>
* </individuals>
* </placements>
* </detector>
* \endcode
*
* @{
*/
// headers
static std::tuple<int, int> add_individuals(Detector& desc, Assembly &env, xml::Collection_t &plm, SensitiveDetector &sens, int id);
static std::tuple<int, int> add_array(Detector& desc, Assembly &env, xml::Collection_t &plm, SensitiveDetector &sens, int id);
static std::tuple<int, int> add_disk(Detector& desc, Assembly &env, xml::Collection_t &plm, SensitiveDetector &sens, int id);
static std::tuple<int, int> 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)
{
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));
std::map<int, int> sectorModuleNumbers;
auto addModuleNumbers = [&sectorModuleNumbers] (int sector, int nmod) {
auto it = sectorModuleNumbers.find(sector);
if (it != sectorModuleNumbers.end()) {
it->second += nmod;
} else {
sectorModuleNumbers[sector] = nmod;
}
};
int sector_id = 1;
for (xml::Collection_t mod(plm, _Unicode(individuals)); mod; ++mod) {
auto [sector, nmod] = add_individuals(desc, assembly, mod, sens, sector_id++);
addModuleNumbers(sector, nmod);
}
for (xml::Collection_t arr(plm, _Unicode(array)); arr; ++arr) {
auto [sector, nmod] = add_array(desc, assembly, arr, sens, sector_id++);
addModuleNumbers(sector, nmod);
}
for (xml::Collection_t disk(plm, _Unicode(disk)); disk; ++disk) {
auto [sector, nmod] = add_disk(desc, assembly, disk, sens, sector_id++);
addModuleNumbers(sector, nmod);
}
for (xml::Collection_t lines(plm, _Unicode(lines)); lines; ++lines) {
auto [sector, nmod] = add_lines(desc, assembly, lines, sens, sector_id++);
addModuleNumbers(sector, nmod);
}
for (auto [sector, nmods] : sectorModuleNumbers) {
desc.add(Constant(Form((detName + "_NModules_Sector%d").c_str(), sector), std::to_string(nmods)));
}
// 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);
envPV.addPhysVolID("system", detID);
det.setPlacement(envPV);
return det;
}
// helper function to build module with or w/o wrapper
std::tuple<Volume, Position> build_module(Detector &desc, xml::Collection_t &plm, SensitiveDetector &sens)
{
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));
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 std::make_tuple(modVol, Position{sx, sy, sz});
// build wrapper
} else {
auto wrp = plm.child(_Unicode(wrapper));
auto thickness = wrp.attr<double>(_Unicode(thickness));
if (thickness < 1e-12*mm) {
return std::make_tuple(modVol, Position{sx, sy, sz});
}
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))));
return std::make_tuple(wrpVol, Position{sx + thickness, sy + thickness, sz});
}
}
// place modules, id must be provided
static std::tuple<int, int> add_individuals(Detector& desc, Assembly &env, xml::Collection_t &plm, SensitiveDetector &sens, int sid)
{
auto [modVol, modSize] = build_module(desc, plm, sens);
int sector_id = dd4hep::getAttrOrDefault<int>(plm, _Unicode(sector), sid);
int nmodules = 0;
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);
nmodules ++;
}
return {sector_id, nmodules};
}
// place array of modules
static std::tuple<int, int> add_array(Detector& desc, Assembly &env, xml::Collection_t &plm, SensitiveDetector &sens, int sid)
{
auto [modVol, modSize] = build_module(desc, plm, sens);
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));
int nmodules = 0;
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);
nmodules ++;
}
}
return {sector_id, nmodules};
}
// place disk of modules
static std::tuple<int, int> add_disk(Detector& desc, Assembly &env, xml::Collection_t &plm, SensitiveDetector &sens, int sid)
{
auto [modVol, modSize] = build_module(desc, plm, sens);
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 = athena::geo::fillRectangles({0., 0.}, modSize.x(), modSize.y(), 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;
for (auto &p : points) {
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);
modPV.addPhysVolID("sector", sector_id).addPhysVolID("module", id_begin + mid++);
}
return {sector_id, mid};
}
// place lines of modules (anchor point is the 0th module of this line)
static std::tuple<int, int> add_lines(Detector& desc, Assembly &env, xml::Collection_t &plm, SensitiveDetector &sens, int sid)
{
auto [modVol, modSize] = build_module(desc, plm, sens);
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);
bool mirrorz = dd4hep::getAttrOrDefault<bool>(plm, _Unicode(mirrorz), false);
// line placement
int mid = 0;
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.));
// determine axis
std::string axis = dd4hep::getAttrOrDefault(pl, _Unicode(axis), "x");
std::transform(axis.begin(), axis.end(), axis.begin(), [](char c) { return std::tolower(c); });
if ((axis != "x") && (axis != "y") && (axis != "z")) {
std::cerr << "HomogeneousCalorimeter Error: cannot determine axis of line " << axis
<< ", abort placement of this line." << std::endl;
continue;
}
int begin = dd4hep::getAttrOrDefault<int>(pl, _Unicode(begin), 0);
int nmods = pl.attr<int>(_Unicode(nmods));
std::vector<Position> trans;
for (int i = 0; i < nmods; ++i) {
Position tran{ (axis == "x") ? pos.x() + (begin + i)*modSize.x() : pos.x(),
(axis == "y") ? pos.y() + (begin + i)*modSize.y() : pos.y(),
(axis == "z") ? pos.z() + (begin + i)*modSize.z() : pos.z() };
trans.push_back(tran);
}
// mirror placement
if (mirrorx) {
int curr_size = trans.size();
for (size_t i = 0; i < curr_size; ++i) {
trans.push_back(Position{-trans[i].x(), trans[i].y(), trans[i].z()});
}
}
if (mirrory) {
int curr_size = trans.size();
for (size_t i = 0; i < curr_size; ++i) {
trans.push_back(Position{trans[i].x(), -trans[i].y(), trans[i].z()});
}
}
if (mirrorz) {
int curr_size = trans.size();
for (size_t i = 0; i < curr_size; ++i) {
trans.push_back(Position{trans[i].x(), trans[i].y(), -trans[i].z()});
}
}
// place volume
for (auto &p : trans) {
Transform3D tr = RotationZYX(rot.z(), rot.y(), rot.x())
* Translation3D(p.x(), p.y(), p.z());
auto modPV = env.placeVolume(modVol, tr);
modPV.addPhysVolID("sector", sector_id).addPhysVolID("module", id_begin + mid++);
}
}
return {sector_id, mid};
}
//@}
DECLARE_DETELEMENT(HomogeneousCalorimeter, create_detector)
src/HybridCalorimeter_geo.cpp 0 → 100644
View file @ 1008a184
//==========================================================================
// 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("SciGlass");
auto crystal_material = desc.material("PbWO4");
auto air_material = desc.material("Air");
double ROut = desc.constantAsDouble("EcalEndcapN_rmax");
double RIn_el = desc.constantAsDouble("EcalEndcapN_rmin1");
double ionCutout_dphi = desc.constantAsDouble("EcalEndcapNIonCutout_dphi");
double RIn = desc.constantAsDouble("EcalEndcapN_rmin2");
double SizeZ = desc.constantAsDouble("EcalEndcapN_thickness");
double thickness = desc.constantAsDouble("EcalEndcapN_thickness");
double trans_radius = desc.constantAsDouble("EcalEndcapNCrystal_rmax");
double Glass_z0 = 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_z0 = desc.constantAsDouble("CrystalModule_z0");
// RIn and ROut will define outer tube embedding the calorimeter
// centers_rmin/out define the maximum radius of module centers
// so that modules are not overlapping with mother tube volume
const double glassHypotenuse = std::hypot(glass_distance, glass_distance)/2;
const double crystalHypotenuse = glassHypotenuse/2;
// Offset these values a bit so we don't miss edge-blocks
const double glassCenters_rmax = ROut - glassHypotenuse + 1 * mm;
const double crystalCenters_rmin = RIn + crystalHypotenuse - 1 * mm;
// Also limits of the inner crystal blocks fill
const double cutout_tan = tan(ionCutout_dphi/2);
const double cutout_rmin = RIn_el + crystalHypotenuse - 1 * mm;
// Offset to align the modules at the zmin of the endcap,
const double Crystal_offset = -0.5 * (Crystal_thickness - thickness);
const double Glass_offset = -0.5*(Glass_thickness - thickness);
// envelope
// consists of an glass tube of the full thickness, and a crystal inner tube
// for the crystal that allows us to get closet to the beampipe without
// overlaps, and a partial electron tube that allows us to get closer to the
// electron beampipe in the region where there is no ion beampipe
Tube glass_tube(min(RIn + glassHypotenuse*2, trans_radius), ROut, SizeZ / 2.0, 0., 360.0 * deg);
Tube crystal_tube(RIn, min(RIn + glassHypotenuse*2, trans_radius), Crystal_thickness/ 2.0, 0., 360.0 * deg);
Tube electron_tube(RIn_el, RIn, Crystal_thickness / 2., ionCutout_dphi / 2., 360.0 * deg - ionCutout_dphi / 2);
UnionSolid outer_envelope(glass_tube, crystal_tube, Position(0,0,Crystal_offset));
UnionSolid envelope(outer_envelope, electron_tube, Position(0,0,Crystal_offset));
Volume ecal_vol("negative_ecal", envelope, air_material);
ecal_vol.setVisAttributes(desc.visAttributes("HybridEcalOuterVis"));
// TODO why 1cm and not something else?
double Glass_OuterR = ROut - 1 * cm ;
double Glass_InnerR = trans_radius;
// 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"));
glass_module.setSensitiveDetector(sens);
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"));
crystal_module.setSensitiveDetector(sens);
// GLASS
double diameter = 2 * Glass_OuterR;
// Can we fit an even or odd amount of glass blocks within our rmax?
// This determines the transition points between crystal and glass as we need the
// outer crystal arangement to be in multiples of 2 (aligned with glass)
const int towersInRow = std::ceil((diameter + Glass_Gap) / (Glass_Width + Glass_Gap));
const bool align_even = (towersInRow % 2);
// 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;
int glass_module_index = 0;
int cryst_module_index = 0;
for(int rowIndex=0; rowIndex < towersInRow; rowIndex++) {
for(int colIndex=0; colIndex < towersInRow; colIndex++) {
const double glass_x = leftTowerPos + colIndex * (Glass_Width + Glass_Gap);
const double glass_y = topTowerPos + rowIndex * (Glass_Width + Glass_Gap);
const double r = std::hypot(glass_x, glass_y);
// crystal if within the transition radius (as defined by the equivalent glass
// block)
if (r < trans_radius) {
for (const auto dx : {-1, 1}) {
for (const auto& dy : {-1, 1}) {
const double crystal_x = glass_x + dx * crystal_distance / 2;
const double crystal_y = glass_y + dy * crystal_distance / 2;
const double crystal_r = std::hypot(crystal_x, crystal_y);
// check if crystal in the main crystal ring?
const bool mainRing = (crystal_r > crystalCenters_rmin);
const bool innerRing = !mainRing && crystal_r > cutout_rmin;
const bool ionCutout = crystal_x > 0 && fabs(crystal_y / crystal_x) < cutout_tan;
if (mainRing || (innerRing && !ionCutout)) {
auto placement =
ecal_vol.placeVolume(crystal_module, Position(crystal_x, crystal_y, Crystal_z0 + Crystal_offset));
placement.addPhysVolID("sector", 1);
placement.addPhysVolID("module", cryst_module_index++);
}
}
}
// Glass block if within the rmax
} else if (r < glassCenters_rmax) {
// glass module
auto placement = ecal_vol.placeVolume(glass_module, Position(glass_x, glass_y, Glass_z0 + Glass_offset));
placement.addPhysVolID("sector", 2);
placement.addPhysVolID("module", glass_module_index++);
}
}
}
desc.add(Constant("EcalEndcapN_NModules_Sector1", std::to_string(cryst_module_index)));
desc.add(Constant("EcalEndcapN_NModules_Sector2", std::to_string(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)
src/MRich_geo.cpp 0 → 100644
View file @ 1008a184
//
// Author : Whit Armstrong (warmstrong@anl.gov)
//
#include <XML/Helper.h>
#include "TMath.h"
#include "TString.h"
#include "DDRec/Surface.h"
#include "DDRec/DetectorData.h"
#include "DD4hep/OpticalSurfaces.h"
#include "DD4hep/DetFactoryHelper.h"
#include "DD4hep/Printout.h"
#include "GeometryHelpers.h"
#include "Math/Vector3D.h"
#include "Math/AxisAngle.h"
#include "Math/VectorUtil.h"
using namespace std;
using namespace dd4hep;
using namespace dd4hep::rec;
using Placements = vector<PlacedVolume>;
static Ref_t createDetector(Detector& description, xml::Handle_t e, SensitiveDetector sens){
xml_det_t x_det = e;
Material air = description.material("AirOptical");
Material vacuum = description.vacuum();
string det_name = x_det.nameStr();
DetElement sdet(det_name, x_det.id());
Assembly assembly(det_name);
sens.setType("tracker");
OpticalSurfaceManager surfMgr = description.surfaceManager();
bool projective = getAttrOrDefault(x_det, _Unicode(projective), false);
bool reflect = x_det.reflect(true);
PlacedVolume pv;
map<string, Volume> modules;
map<string, Placements> sensitives;
map<string, Volume> module_assemblies;
std::map<std::string,DetElement> module_assembly_delements;
int n_sensor = 1;
// dimensions
xml::Component dims = x_det.dimensions();
auto rmin = dims.rmin();
auto rmax = dims.rmax();
auto length = dims.length();
auto zmin = dims.zmin();
auto zpos = zmin + length / 2;
// envelope
Tube envShape(rmin, rmax, length / 2., 0., 2 * M_PI);
Volume envVol("MRICH_Envelope", envShape, air);
envVol.setVisAttributes(description.visAttributes(x_det.visStr()));
if (x_det.hasChild(_Unicode(envelope))) {
xml_comp_t x_envelope = x_det.child(_Unicode(envelope));
double thickness = x_envelope.thickness();
Material material = description.material(x_envelope.materialStr());
Tube envInsideShape(rmin + thickness, rmax - thickness, length / 2. - thickness);
SubtractionSolid envShellShape(envShape, envInsideShape);
Volume envShell("MRICH_Envelope_Inside", envShellShape, material);
envVol.placeVolume(envShell);
}
// expect only one module (for now)
xml_comp_t x_mod = x_det.child(_U(module));
string mod_name = x_mod.nameStr();
double mod_width = getAttrOrDefault(x_mod, _U(width), 130.0 * mm);
double mod_height = getAttrOrDefault(x_mod, _U(height), 130.0 * mm);
double mod_length = getAttrOrDefault(x_mod, _U(length), 130.0 * mm);
// module
Box m_solid(mod_width / 2.0, mod_height / 2.0, mod_length / 2.0);
Volume m_volume(mod_name, m_solid, air);
m_volume.setVisAttributes(description.visAttributes(x_mod.visStr()));
DetElement mod_de( mod_name + std::string("_mod_") + std::to_string(1), 1);
double z_placement = - mod_length / 2.0;
// todo module frame
if (x_mod.hasChild(_Unicode(frame))) {
xml_comp_t x_frame = x_mod.child(_Unicode(frame));
double frame_thickness = getAttrOrDefault(x_frame, _U(thickness), 2.0 * mm);
Box frame_inside(mod_width / 2.0 - frame_thickness, mod_height / 2.0 - frame_thickness, mod_length / 2.0 - frame_thickness);
SubtractionSolid frame_solid(m_solid, frame_inside);
Material frame_mat = description.material(x_frame.materialStr());
Volume frame_vol(mod_name+"_frame", frame_solid, frame_mat);
auto frame_vis = getAttrOrDefault<std::string>(x_frame, _U(vis), std::string("GrayVis"));
frame_vol.setVisAttributes(description.visAttributes(frame_vis));
// update position
z_placement += frame_thickness / 2.0;
// place volume
m_volume.placeVolume(frame_vol);
// update position
z_placement += frame_thickness / 2.0;
}
// aerogel box
if (x_mod.hasChild(_Unicode(aerogel))) {
xml_comp_t x_aerogel = x_mod.child(_Unicode(aerogel));
double aerogel_width = getAttrOrDefault(x_aerogel, _U(width), 130.0 * mm);
double aerogel_length = getAttrOrDefault(x_aerogel, _U(length), 130.0 * mm);
Material aerogel_mat = description.material(x_aerogel.materialStr());
auto aerogel_vis = getAttrOrDefault<std::string>(x_aerogel, _U(vis), std::string("InvisibleWithDaughters"));
xml_comp_t x_aerogel_frame = x_aerogel.child(_Unicode(frame));
double foam_thickness = getAttrOrDefault(x_aerogel_frame, _U(thickness), 2.0 * mm);
Material foam_mat = description.material(x_aerogel_frame.materialStr());
auto foam_vis = getAttrOrDefault<std::string>(x_aerogel_frame, _U(vis), std::string("RedVis"));
// foam frame
Box foam_box(aerogel_width / 2.0 + foam_thickness, aerogel_width / 2.0 + foam_thickness, (aerogel_length + foam_thickness) / 2.0);
Box foam_sub_box(aerogel_width / 2.0, aerogel_width / 2.0, (aerogel_length + foam_thickness) / 2.0);
SubtractionSolid foam_frame_solid(foam_box, foam_sub_box, Position(0, 0, foam_thickness));
Volume foam_vol(mod_name+"_aerogel_frame", foam_frame_solid, foam_mat);
foam_vol.setVisAttributes(description.visAttributes(foam_vis));
// aerogel
Box aerogel_box(aerogel_width / 2.0, aerogel_width / 2.0, (aerogel_length) / 2.0);
Volume aerogel_vol(mod_name+"_aerogel", aerogel_box, aerogel_mat);
aerogel_vol.setVisAttributes(description.visAttributes(aerogel_vis));
// update position
z_placement += (aerogel_length + foam_thickness) / 2.0;
// place foam frame
pv = m_volume.placeVolume(foam_vol,Position(0,0,z_placement));
// place aerogel
z_placement += foam_thickness / 2.0;
pv = m_volume.placeVolume(aerogel_vol,Position(0,0,z_placement));
DetElement aerogel_de(mod_de, mod_name + std::string("_aerogel_de") + std::to_string(1), 1);
aerogel_de.setPlacement(pv);
// update position
z_placement += aerogel_length / 2.0;
// optical surfaces
auto aerogel_surf = surfMgr.opticalSurface(dd4hep::getAttrOrDefault(x_aerogel, _Unicode(surface), "MRICH_AerogelOpticalSurface"));
SkinSurface skin_surf(description, aerogel_de, Form("MRICH_aerogel_skin_surface_%d", 1), aerogel_surf, aerogel_vol);
skin_surf.isValid();
}
// Fresnel Lens
if (x_mod.hasChild(_Unicode(lens))) {
xml_comp_t x_lens = x_mod.child(_Unicode(lens));
// - The lens has a constant groove pitch (delta r) as opposed to fixing the groove height.
// - The lens area outside of the effective diamtere is flat.
// - The grooves are not curved, rather they are polycone shaped, ie a flat approximating the curvature.
auto lens_vis = getAttrOrDefault<std::string>(x_lens, _U(vis), std::string("AnlBlue"));
double groove_pitch = getAttrOrDefault(x_lens, _Unicode(pitch), 0.2 * mm);// 0.5 * mm);
double lens_f = getAttrOrDefault(x_lens, _Unicode(focal_length), 6.0*2.54*cm);
double eff_diameter = getAttrOrDefault(x_lens, _Unicode(effective_diameter), 152.4 * mm);
double lens_width = getAttrOrDefault(x_lens, _Unicode(width), 6.7*2.54*cm);
double center_thickness = getAttrOrDefault(x_lens, _U(thickness), 0.068 * 2.54 * cm);//2.0 * mm);
double n_acrylic = 1.49;
double lens_curvature = 1.0 / (lens_f*(n_acrylic - 1.0)); //confirmed
double full_ring_rmax = std::min(eff_diameter / 2.0, lens_width/2.0);
double N_grooves = std::ceil((full_ring_rmax) / groove_pitch);
double groove_last_rmin = (N_grooves - 1) * groove_pitch;
double groove_last_rmax = N_grooves * groove_pitch;
auto groove_sagitta = [&](double r) { return lens_curvature * std::pow(r, 2) / (1.0 + 1.0); };
double lens_thickness = groove_sagitta(groove_last_rmax) - groove_sagitta(groove_last_rmin) + center_thickness;
Material lens_mat = description.material(x_lens.materialStr());
Box lens_box(lens_width / 2.0, lens_width / 2.0, (center_thickness) / 2.0);
SubtractionSolid flat_lens(lens_box, Tube(0.0, full_ring_rmax, 2 * center_thickness));
Assembly lens_vol(mod_name + "_lens");
Volume flatpart_lens_vol( "flatpart_lens", flat_lens, lens_mat);
lens_vol.placeVolume(flatpart_lens_vol);
int i_groove = 0;
double groove_rmax = groove_pitch;
double groove_rmin = 0;
while ( groove_rmax <= full_ring_rmax ) {
double dZ = groove_sagitta(groove_rmax) - groove_sagitta(groove_rmin);
Polycone groove_solid(0, 2.0 * M_PI,
{groove_rmin, groove_rmin, groove_rmin},
{groove_rmax, groove_rmax, groove_rmin},
{-lens_thickness/2.0, lens_thickness/2.0-dZ, lens_thickness/2.0});
Volume lens_groove_vol("lens_groove_" + std::to_string(i_groove), groove_solid, lens_mat);
lens_vol.placeVolume(lens_groove_vol);
i_groove++;
groove_rmin = (i_groove )*groove_pitch;
groove_rmax = (i_groove+1)*groove_pitch;
}
lens_vol.setVisAttributes(description.visAttributes(lens_vis));
// update position
z_placement += lens_thickness/2.0;
// place volume
pv = m_volume.placeVolume(lens_vol,Position(0,0,z_placement));
DetElement lens_de(mod_de, mod_name + std::string("_lens_de") + std::to_string(1), 1);
lens_de.setPlacement(pv);
// update position
z_placement += lens_thickness/2.0;
// optical surfaces
auto lens_surf = surfMgr.opticalSurface(dd4hep::getAttrOrDefault(x_lens, _Unicode(surface), "MRICH_LensOpticalSurface"));
SkinSurface skin_surf(description, lens_de, Form("MRichFresnelLens_skin_surface_%d", 1), lens_surf, lens_vol);
skin_surf.isValid();
}
// mirror
if (x_mod.hasChild(_Unicode(space))) {
xml_comp_t x_space = x_mod.child(_Unicode(space));
z_placement += getAttrOrDefault(x_space, _U(thickness), 0.0 * mm);
}
// mirror
if (x_mod.hasChild(_Unicode(mirror))) {
xml_comp_t x_mirror = x_mod.child(_Unicode(mirror));
auto mirror_vis = getAttrOrDefault<std::string>(x_mirror, _U(vis), std::string("AnlGray"));
double mirror_x1 = getAttrOrDefault(x_mirror, _U(x1), 100.0 * mm);
double mirror_x2 = getAttrOrDefault(x_mirror, _U(x2), 80.0 * mm);
double mirror_length = getAttrOrDefault(x_mirror, _U(length), 130.0 * mm);
double mirror_thickness = getAttrOrDefault(x_mirror, _U(thickness), 2.0 * mm);
double outer_x1 = (mirror_x1+mirror_thickness)/2.0;
double outer_x2 = (mirror_x2+mirror_thickness)/2.0;
Trd2 outer_mirror_trd(outer_x1, outer_x2, outer_x1, outer_x2, mirror_length/2.0);
Trd2 inner_mirror_trd(mirror_x1 / 2.0, mirror_x2 / 2.0, mirror_x1 / 2.0,mirror_x2 / 2.0, mirror_length/2.0+0.1*mm);
SubtractionSolid mirror_solid(outer_mirror_trd, inner_mirror_trd);
Material mirror_mat = description.material(x_mirror.materialStr());
Volume mirror_vol(mod_name+"_mirror", mirror_solid, mirror_mat);
// update position
z_placement += mirror_length/2.0;
// place volume
pv = m_volume.placeVolume(mirror_vol,Position(0,0,z_placement));
DetElement mirror_de(mod_de, mod_name + std::string("_mirror_de") + std::to_string(1), 1);
mirror_de.setPlacement(pv);
// update position
z_placement += mirror_length/2.0;
// optical surfaces
auto mirror_surf = surfMgr.opticalSurface(dd4hep::getAttrOrDefault(x_mirror, _Unicode(surface), "MRICH_MirrorOpticalSurface"));
SkinSurface skin_surf(description, mirror_de, Form("MRICH_mirror_skin_surface_%d", 1), mirror_surf, mirror_vol);
skin_surf.isValid();
}
// photon detector
if (x_mod.hasChild(_Unicode(photodet))) {
xml_comp_t x_photodet = x_mod.child(_Unicode(photodet));
auto photodet_vis = getAttrOrDefault<std::string>(x_photodet, _U(vis), std::string("AnlRed"));
double photodet_width = getAttrOrDefault(x_photodet, _U(width), 130.0 * mm);
double photodet_thickness = getAttrOrDefault(x_photodet, _U(thickness), 2.0 * mm);
Material photodet_mat = description.material(x_photodet.materialStr());
Box window_box(photodet_width/2.0,photodet_width/2.0,photodet_thickness/2.0);
Volume window_vol(mod_name+"_window", window_box, photodet_mat);
// update position
z_placement += photodet_thickness/2.0;
// place volume
pv = m_volume.placeVolume(window_vol,Position(0,0,z_placement));
DetElement comp_de(mod_de, mod_name + std::string("_sensor_de_") + std::to_string(1) , 1);
comp_de.setPlacement(pv);
// update position
z_placement += photodet_thickness/2.0;
// sensitive
pv.addPhysVolID("sensor", n_sensor);
window_vol.setSensitiveDetector(sens);
sensitives[mod_name].push_back(pv);
++n_sensor;
// sensor
xml_comp_t x_sensor = x_photodet.child(_Unicode(sensor));
double sensor_thickness = getAttrOrDefault(x_sensor, _U(thickness), 2.0 * mm);
Material sensor_mat = description.material(x_sensor.materialStr());
int sensor_nx = getAttrOrDefault(x_sensor, _Unicode(nx), 2);
int sensor_ny = getAttrOrDefault(x_sensor, _Unicode(ny), 2);
// layers
int i_layer = 1;
for (xml_coll_t li(x_photodet, _Unicode(layer)); li; ++li) {
xml_comp_t x_layer = li;
Material layer_mat = description.material(x_layer.materialStr());
double layer_thickness = x_layer.thickness();
Box layer_box(photodet_width/2.0,photodet_width/2.0,layer_thickness/2.0);
Volume layer_vol(mod_name + "_layer_" + std::to_string(i_layer), layer_box, layer_mat);
// update position
z_placement += layer_thickness / 2.0;
// place volume
pv = m_volume.placeVolume(layer_vol,Position(0,0,z_placement));
DetElement layer_de(mod_de, mod_name + std::string("_layer_de_") + std::to_string(i_layer), 1);
layer_de.setPlacement(pv);
// update position
z_placement += layer_thickness / 2.0;
i_layer++;
}
}
//for (size_t ic = 0; ic < sensVols.size(); ++ic) {
// PlacedVolume sens_pv = sensVols[ic];
// DetElement comp_de(mod_de, std::string("de_") + sens_pv.volume().name(), ic + 1);
// comp_de.setPlacement(sens_pv);
// // Acts::ActsExtension* sensorExtension = new Acts::ActsExtension();
// //// sensorExtension->addType("sensor", "detector");
// // comp_de.addExtension<Acts::ActsExtension>(sensorExtension);
// //// comp_de.setAttributes(description, sens_pv.volume(), x_layer.regionStr(),
// //// x_layer.limitsStr(),
// //// xml_det_t(xmleles[m_nam]).visStr());
//}
//DetElement window_de(sdet, mod_name + std::string("_window_de") + std::to_string(1), 1);
//window_de.setPlacement(pv);
modules[mod_name] = m_volume;
module_assembly_delements[mod_name] = mod_de;
// end module
// place modules in the sectors (disk)
auto points = athena::geo::fillSquares({0., 0.}, mod_width, rmin, rmax);
// mod_name = ...
Placements& sensVols = sensitives[mod_name];
auto mod_v = modules[mod_name];
// determine module direction, always facing z = 0
double roty = dims.zmin() < 0. ? -M_PI : 0 ;
// read module positions
std::vector<std::tuple<double,double,double>> positions;
for (xml_coll_t x_positions_i(x_det, _Unicode(positions)); x_positions_i; ++x_positions_i) {
xml_comp_t x_positions = x_positions_i;
for (xml_coll_t x_position_i(x_positions, _U(position)); x_position_i; ++x_position_i) {
xml_comp_t x_position = x_position_i;
positions.push_back(
std::make_tuple(x_positions.scale() * x_position.x() * mm,
x_positions.scale() * x_position.y() * mm,
-x_positions.z0()));
}
}
// if no positions, then autoplacement
if (positions.empty()) {
for (double x = mod_width / 2.0; x < rmax - mod_width / 2.0; x += mod_width) {
for (double y = mod_width / 2.0; y < rmax - mod_width / 2.0; y += mod_width) {
if (pow(x + mod_width / 2.0,2) + pow(y + mod_width / 2.0,2) > rmax*rmax) continue;
if (pow(x - mod_width / 2.0,2) + pow(y - mod_width / 2.0,2) < rmin*rmin) continue;
positions.push_back(std::make_tuple(x, y, 0));
}
}
}
// place modules
int i_mod = 1; // starts at 1
for (auto& p: positions) {
// get positions in one quadrant
double x = std::get<0>(p);
double y = std::get<1>(p);
double z0 = std::get<2>(p);
// and place in all quadrants (intentional shadowing)
for (auto& p: decltype(positions){{x,y,z0}, {y,-x,z0}, {-x,-y,z0}, {-y,x,z0}}) {
// get positions (intentional shadowing)
double x = std::get<0>(p);
double y = std::get<1>(p);
double z0 = std::get<2>(p);
Transform3D tr;
if(projective) {
double rotAngX = atan(y/z0);
double rotAngY = -1.*atan(x/z0);
tr = Translation3D(x, y, 0) * RotationX(rotAngX) * RotationY(rotAngY);
} else {
tr = Translation3D(x, y, 0) * RotationX(0);
}
// mod placement
pv = envVol.placeVolume(mod_v, tr);
pv.addPhysVolID("module", i_mod);
auto mod_det_element = module_assembly_delements[mod_name].clone(mod_name + "__" + std::to_string(i_mod));
mod_det_element.setPlacement(pv);
sdet.add(mod_det_element);
i_mod++;
}
}
// additional layers
if (x_det.hasChild(_Unicode(layer))) {
xml_comp_t x_layer = x_det.child(_Unicode(layer));
double layer_thickness = x_layer.thickness();
Material layer_mat = description.material(x_layer.materialStr());
Tube frameShape(rmin, rmax, layer_thickness / 2., 0., 2 * M_PI);
Volume frameVol("MRICH_Frame", frameShape, layer_mat);
pv = envVol.placeVolume(frameVol, Position(0, 0, (length - layer_thickness) / 2.0));
}
// place envelope
Volume motherVol = description.pickMotherVolume(sdet);
if (reflect) {
pv = motherVol.placeVolume(envVol, Transform3D(RotationZYX(0, M_PI, 0), Position(0, 0, -zpos)));
} else {
pv = motherVol.placeVolume(envVol, Transform3D(RotationZYX(0, 0, 0), Position(0, 0, +zpos)));
}
pv.addPhysVolID("system", x_det.id());
sdet.setPlacement(pv);
return sdet;
}
//void addModules(Volume &mother, xml::DetElement &detElem, Detector &description, SensitiveDetector &sens)
//{
// xml::Component dims = detElem.dimensions();
// xml::Component mods = detElem.child(_Unicode(modules));
//
// auto rmin = dims.rmin();
// auto rmax = dims.rmax();
//
// auto mThick = mods.attr<double>(_Unicode(thickness));
// auto mWidth = mods.attr<double>(_Unicode(width));
// auto mGap = mods.attr<double>(_Unicode(gap));
//
// auto modMat = description.material(mods.materialStr());
// auto gasMat = description.material("AirOptical");
//
// // single module
// Box mShape(mWidth/2., mWidth/2., mThick/2. - 0.1*mm);
// Volume mVol("ce_MRICH_mod_Solid", mShape, modMat);
//
// // a thin gas layer to detect optical photons
// Box modShape(mWidth/2., mWidth/2., mThick/2.);
// Volume modVol("ce_MRICH_mod_Solid_v", modShape, gasMat);
// // thin gas layer is on top (+z) of the material
// modVol.placeVolume(mVol, Position(0., 0., -0.1*mm));
//
// modVol.setVisAttributes(description.visAttributes(mods.visStr()));
// sens.setType("tracker");
// modVol.setSensitiveDetector(sens);
//
// // place modules in the sectors (disk)
// auto points = ref::utils::fillSquares({0., 0.}, mWidth + mGap, rmin - mGap, rmax + mGap);
//
// // determine module direction, always facing z = 0
// double roty = dims.z() > 0. ? M_PI/2. : -M_PI/2.;
// int imod = 1;
// for (auto &p : points) {
// // operations are inversely ordered
// Transform3D tr = Translation3D(p.x(), p.y(), 0.) // move to position
// * RotationY(roty); // facing z = 0.
// auto modPV = mother.placeVolume(modVol, tr);
// modPV.addPhysVolID("sector", 1).addPhysVolID("module", imod ++);
// }
//}
// clang-format off
DECLARE_DETELEMENT(athena_MRICH, createDetector)
src/PFRICH_geo.cpp 0 → 100644
View file @ 1008a184
//----------------------------------
// pfRICH: Proximity Focusing RICH
// Author: C. Dilks
//----------------------------------
#include "DD4hep/DetFactoryHelper.h"
#include "DD4hep/OpticalSurfaces.h"
#include "DD4hep/Printout.h"
#include "DDRec/DetectorData.h"
#include "DDRec/Surface.h"
#include "GeometryHelpers.h"
#include "Math/Point2D.h"
#include "TMath.h"
#include "TString.h"
#include <XML/Helper.h>
using namespace dd4hep;
using namespace dd4hep::rec;
// create the detector
static Ref_t createDetector(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);
xml::Component dims = detElem.dimensions();
OpticalSurfaceManager surfMgr = desc.surfaceManager();
// attributes -----------------------------------------------------------
// - vessel
double vesselLength = dims.attr<double>(_Unicode(length));
double vesselZmin = dims.attr<double>(_Unicode(zmin));
double vesselZmax = dims.attr<double>(_Unicode(zmax));
double vesselRmin0 = dims.attr<double>(_Unicode(rmin0));
double vesselRmin1 = dims.attr<double>(_Unicode(rmin1));
double vesselRmax0 = dims.attr<double>(_Unicode(rmax0));
double vesselRmax1 = dims.attr<double>(_Unicode(rmax1));
double wallThickness = dims.attr<double>(_Unicode(wall_thickness));
double windowThickness = dims.attr<double>(_Unicode(window_thickness));
auto vesselMat = desc.material(detElem.attr<std::string>(_Unicode(material)));
auto gasvolMat = desc.material(detElem.attr<std::string>(_Unicode(gas)));
auto vesselVis = desc.visAttributes(detElem.attr<std::string>(_Unicode(vis_vessel)));
auto gasvolVis = desc.visAttributes(detElem.attr<std::string>(_Unicode(vis_gas)));
// - radiator (applies to aerogel and filter)
auto radiatorElem = detElem.child(_Unicode(radiator));
double radiatorRmin = radiatorElem.attr<double>(_Unicode(rmin));
double radiatorRmax = radiatorElem.attr<double>(_Unicode(rmax));
double radiatorPhiw = radiatorElem.attr<double>(_Unicode(phiw));
double radiatorPitch = radiatorElem.attr<double>(_Unicode(pitch));
double radiatorFrontplane = radiatorElem.attr<double>(_Unicode(frontplane));
// - aerogel
auto aerogelElem = radiatorElem.child(_Unicode(aerogel));
auto aerogelMat = desc.material(aerogelElem.attr<std::string>(_Unicode(material)));
auto aerogelVis = desc.visAttributes(aerogelElem.attr<std::string>(_Unicode(vis)));
double aerogelThickness = aerogelElem.attr<double>(_Unicode(thickness));
// - filter
auto filterElem = radiatorElem.child(_Unicode(filter));
auto filterMat = desc.material(filterElem.attr<std::string>(_Unicode(material)));
auto filterVis = desc.visAttributes(filterElem.attr<std::string>(_Unicode(vis)));
double filterThickness = filterElem.attr<double>(_Unicode(thickness));
// - sensor module
auto sensorElem = detElem.child(_Unicode(sensors)).child(_Unicode(module));
auto sensorMat = desc.material(sensorElem.attr<std::string>(_Unicode(material)));
auto sensorVis = desc.visAttributes(sensorElem.attr<std::string>(_Unicode(vis)));
auto sensorSurf = surfMgr.opticalSurface(sensorElem.attr<std::string>(_Unicode(surface)));
double sensorSide = sensorElem.attr<double>(_Unicode(side));
double sensorGap = sensorElem.attr<double>(_Unicode(gap));
double sensorThickness = sensorElem.attr<double>(_Unicode(thickness));
// - sensor plane
auto sensorPlaneElem = detElem.child(_Unicode(sensors)).child(_Unicode(plane));
double sensorPlaneDist = sensorPlaneElem.attr<double>(_Unicode(sensordist));
double sensorPlaneRmin = sensorPlaneElem.attr<double>(_Unicode(rmin));
double sensorPlaneRmax = sensorPlaneElem.attr<double>(_Unicode(rmax));
// - debugging switches
int debug_optics_mode = detElem.attr<int>(_Unicode(debug_optics));
// if debugging optics, override some settings
bool debug_optics = debug_optics_mode > 0;
if (debug_optics) {
printout(WARNING, "PFRICH_geo", "DEBUGGING PFRICH OPTICS");
switch (debug_optics_mode) {
case 1:
vesselMat = aerogelMat = filterMat = sensorMat = gasvolMat = desc.material("VacuumOptical");
break;
case 2:
vesselMat = aerogelMat = filterMat = sensorMat = desc.material("VacuumOptical");
break;
default:
printout(FATAL, "PFRICH_geo", "UNKNOWN debug_optics_mode");
return det;
};
aerogelVis = sensorVis;
gasvolVis = vesselVis = desc.invisible();
};
// BUILD VESSEL //////////////////////////////////////
/* - `vessel`: aluminum enclosure, the mother volume of the pfRICH
* - `gasvol`: gas volume, which fills `vessel`; all other volumes defined below
* are children of `gasvol`
*/
// tank solids
double boreDelta = vesselRmin1 - vesselRmin0;
Cone vesselTank(vesselLength / 2.0, vesselRmin1, vesselRmax1, vesselRmin0, vesselRmax0);
Cone gasvolTank(vesselLength / 2.0 - windowThickness, vesselRmin1 + wallThickness, vesselRmax1 - wallThickness,
vesselRmin0 + wallThickness, vesselRmax0 - wallThickness);
// extra solids for `debug_optics` only
Box vesselBox(1001, 1001, 1001);
Box gasvolBox(1000, 1000, 1000);
// choose vessel and gasvol solids (depending on `debug_optics_mode` (0=disabled))
Solid vesselSolid, gasvolSolid;
switch (debug_optics_mode) {
case 0:
vesselSolid = vesselTank;
gasvolSolid = gasvolTank;
break; // `!debug_optics`
case 1:
vesselSolid = vesselBox;
gasvolSolid = gasvolBox;
break;
case 2:
vesselSolid = vesselBox;
gasvolSolid = gasvolTank;
break;
};
// volumes
Volume vesselVol(detName, vesselSolid, vesselMat);
Volume gasvolVol(detName + "_gas", gasvolSolid, gasvolMat);
vesselVol.setVisAttributes(vesselVis);
gasvolVol.setVisAttributes(gasvolVis);
// reference positions
// - the vessel is created such that the center of the cylindrical tank volume
// coincides with the origin; this is called the "origin position" of the vessel
// - when the vessel (and its children volumes) is placed, it is translated in
// the z-direction to be in the proper ATHENA-integration location
// - these reference positions are for the frontplane and backplane of the vessel,
// with respect to the vessel origin position
auto originFront = Position(0., 0., vesselLength / 2.0);
auto originBack = Position(0., 0., -vesselLength / 2.0);
// sensitive detector type
sens.setType("tracker");
// BUILD RADIATOR //////////////////////////////////////
// attributes
double airGap = 0.01*mm; // air gap between aerogel and filter (FIXME? actually it's currently a gas gap)
// solid and volume: create aerogel and filter
Cone aerogelSolid(
aerogelThickness/2,
radiatorRmin + boreDelta * aerogelThickness / vesselLength, /* at backplane */
radiatorRmax,
radiatorRmin, /* at frontplane */
radiatorRmax
);
Cone filterSolid(
filterThickness/2,
radiatorRmin + boreDelta * (aerogelThickness + airGap + filterThickness) / vesselLength, /* at backplane */
radiatorRmax,
radiatorRmin + boreDelta * (aerogelThickness + airGap) / vesselLength, /* at frontplane */
radiatorRmax
);
Volume aerogelVol(detName + "_aerogel", aerogelSolid, aerogelMat);
Volume filterVol(detName + "_filter", filterSolid, filterMat);
aerogelVol.setVisAttributes(aerogelVis);
filterVol.setVisAttributes(filterVis);
// aerogel placement and surface properties
// TODO [low-priority]: define skin properties for aerogel and filter
auto radiatorPos = Position(0., 0., radiatorFrontplane - 0.5 * aerogelThickness) + originFront;
auto aerogelPV = gasvolVol.placeVolume(
aerogelVol,
Translation3D(radiatorPos.x(), radiatorPos.y(), radiatorPos.z()) // re-center to originFront
* RotationY(radiatorPitch) // change polar angle to specified pitch // (FIXME: probably broken, currently not in use)
);
DetElement aerogelDE(det, "aerogel_de", 0);
aerogelDE.setPlacement(aerogelPV);
// SkinSurface aerogelSkin(desc, aerogelDE, "mirror_optical_surface", aerogelSurf, aerogelVol);
// aerogelSkin.isValid();
// filter placement and surface properties
if (!debug_optics) {
auto filterPV = gasvolVol.placeVolume(
filterVol,
Translation3D(0., 0., -airGap) // add an airgap (FIXME: actually a gas gap)
* Translation3D(radiatorPos.x(), radiatorPos.y(), radiatorPos.z()) // re-center to originFront
* RotationY(radiatorPitch) // change polar angle
* Translation3D(0., 0., -(aerogelThickness + filterThickness) / 2.) // move to aerogel backplane
);
DetElement filterDE(det, "filter_de", 0);
filterDE.setPlacement(filterPV);
// SkinSurface filterSkin(desc, filterDE, "mirror_optical_surface", filterSurf, filterVol);
// filterSkin.isValid();
};
// BUILD SENSORS ///////////////////////
// solid and volume: single sensor module
Box sensorSolid(sensorSide / 2., sensorSide / 2., sensorThickness / 2.);
Volume sensorVol(detName + "_sensor", sensorSolid, sensorMat);
sensorVol.setVisAttributes(sensorVis);
// sensitivity
if (!debug_optics)
sensorVol.setSensitiveDetector(sens);
// sensor plane positioning: we want `sensorPlaneDist` to be the distance between the
// aerogel backplane (i.e., aerogel/filter boundary) and the sensor active surface (e.g, photocathode)
double sensorZpos = radiatorFrontplane - aerogelThickness - sensorPlaneDist - 0.5 * sensorThickness;
auto sensorPlanePos = Position(0., 0., sensorZpos) + originFront; // reference position
// miscellaneous
int imod = 0; // module number
double tBoxMax = vesselRmax1; // sensors will be tiled in tBox, within annular limits
// SENSOR MODULE LOOP ------------------------
/* cartesian tiling loop
* - start at (x=0,y=0), to center the grid
* - loop over positive-x positions; for each, place the corresponding negative-x sensor too
* - nested similar loop over y positions
*/
double sx, sy;
for (double usx = 0; usx <= tBoxMax; usx += sensorSide + sensorGap) {
for (int sgnx = 1; sgnx >= (usx > 0 ? -1 : 1); sgnx -= 2) {
for (double usy = 0; usy <= tBoxMax; usy += sensorSide + sensorGap) {
for (int sgny = 1; sgny >= (usy > 0 ? -1 : 1); sgny -= 2) {
// sensor (x,y) center
sx = sgnx * usx;
sy = sgny * usy;
// annular cut
if (std::hypot(sx, sy) < sensorPlaneRmin || std::hypot(sx, sy) > sensorPlaneRmax)
continue;
// placement (note: transformations are in reverse order)
auto sensorPV = gasvolVol.placeVolume(
sensorVol, Transform3D(Translation3D(sensorPlanePos.x(), sensorPlanePos.y(),
sensorPlanePos.z()) // move to reference position
* Translation3D(sx, sy, 0.) // move to grid position
));
// generate LUT for module number -> sensor position, for readout mapping tests
// printf("%d %f %f\n",imod,sensorPV.position().x(),sensorPV.position().y());
// properties
sensorPV.addPhysVolID("module", imod);
DetElement sensorDE(det, Form("sensor_de_%d", imod), imod);
sensorDE.setPlacement(sensorPV);
if (!debug_optics) {
SkinSurface sensorSkin(desc, sensorDE, "sensor_optical_surface", sensorSurf,
sensorVol); // TODO: 3rd arg needs `imod`?
sensorSkin.isValid();
};
// increment sensor module number
imod++;
};
};
};
};
// END SENSOR MODULE LOOP ------------------------
//
// Add service material if desired
if (detElem.child("sensors").hasChild("services")) {
xml_comp_t x_service = detElem.child("sensors").child(_Unicode(services));
Assembly service_vol("services");
service_vol.setVisAttributes(desc, x_service.visStr());
// Compute service total thickness from components
double total_thickness = 0;
xml_coll_t ci(x_service, _Unicode(component));
for (ci.reset(), total_thickness = 0.0; ci; ++ci) {
total_thickness += xml_comp_t(ci).thickness();
}
int ncomponents = 0;
double thickness_sum = -total_thickness / 2.0;
for (xml_coll_t ci(x_service, _Unicode(component)); ci; ++ci, ncomponents++) {
xml_comp_t x_comp = ci;
double thickness = x_comp.thickness();
Tube c_tube{sensorPlaneRmin, sensorPlaneRmax, thickness/2};
Volume c_vol{_toString(ncomponents, "component%d"), c_tube, desc.material(x_comp.materialStr())};
c_vol.setVisAttributes(desc, x_comp.visStr());
service_vol.placeVolume(c_vol, Position(0, 0, thickness_sum + thickness / 2.0));
thickness_sum += thickness;
}
gasvolVol.placeVolume(service_vol,
Transform3D(Translation3D(sensorPlanePos.x(), sensorPlanePos.y(),
sensorPlanePos.z() - sensorThickness - total_thickness)));
}
// place gas volume
PlacedVolume gasvolPV = vesselVol.placeVolume(gasvolVol, Position(0, 0, 0));
DetElement gasvolDE(det, "gasvol_de", 0);
gasvolDE.setPlacement(gasvolPV);
// place mother volume (vessel)
Volume motherVol = desc.pickMotherVolume(det);
PlacedVolume vesselPV = motherVol.placeVolume(vesselVol, Position(0, 0, vesselZmin) - originFront);
vesselPV.addPhysVolID("system", detID);
det.setPlacement(vesselPV);
return det;
};
// clang-format off
DECLARE_DETELEMENT(athena_PFRICH, createDetector)
src/PolyhedraEndcapCalorimeter2_geo.cpp
View file @ 1008a184
...@@ -11,7 +11,7 @@ ...@@ -11,7 +11,7 @@
// //
//========================================================================== //==========================================================================
// //
// Modified for TOPSiDE detector // Modified for ATHENA detector
// //
//========================================================================== //==========================================================================
#include "DD4hep/DetFactoryHelper.h" #include "DD4hep/DetFactoryHelper.h"
...@@ -32,7 +32,7 @@ static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector s ...@@ -32,7 +32,7 @@ static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector s
int numsides = dim.numsides(); int numsides = dim.numsides();
xml::Component pos = x_det.position(); xml::Component pos = x_det.position();
double rmin = dim.rmin(); double rmin = dim.rmin();
double rmax = dim.rmax() * std::cos(M_PI / numsides); double rmax = dim.rmax();
double zmin = dim.zmin(); double zmin = dim.zmin();
Layering layering(x_det); Layering layering(x_det);
double totalThickness = layering.totalThickness(); double totalThickness = layering.totalThickness();
...@@ -128,4 +128,5 @@ static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector s ...@@ -128,4 +128,5 @@ static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector s
} }
// clang-format off // clang-format off
DECLARE_DETELEMENT(refdet_PolyhedraEndcapCalorimeter2, create_detector) DECLARE_DETELEMENT(athena_PolyhedraEndcapCalorimeter2, create_detector)
DECLARE_DETELEMENT(athena_PolyhedraEndcapCalorimeter, create_detector)
src/ScFiCalorimeter_geo.cpp 0 → 100644
View file @ 1008a184
//==========================================================================
// 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()));
// build module
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});
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))));
}
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);
}
}
// if no fibers we make the module itself sensitive
} else {
modVol.setSensitiveDetector(sens);
}
return std::make_tuple(modVol, Position{sx, sy, sz});
}
DECLARE_DETELEMENT(ScFiCalorimeter, create_detector)
src/ShashlikCalorimeter_geo.cpp 0 → 100644
View file @ 1008a184
//==========================================================================
// Implementation for shashlik calorimeter modules
// it supports disk placements with (rmin, rmax), and (phimin, phimax)
//--------------------------------------------------------------------------
// Author: Chao Peng (ANL)
// Date: 06/22/2021
//==========================================================================
#include "GeometryHelpers.h"
#include "DD4hep/DetFactoryHelper.h"
#include <XML/Layering.h>
#include <XML/Helper.h>
#include <iostream>
#include <algorithm>
#include <tuple>
#include <math.h>
using namespace dd4hep;
static void add_disk_shashlik(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;
}
static Ref_t create_detector(Detector& desc, xml::Handle_t handle, SensitiveDetector sens)
{
static const std::string func = "ShashlikCalorimeter";
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(disk)); mod; ++mod) {
add_disk_shashlik(desc, assembly, mod, 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);
envPV.addPhysVolID("system", detID);
det.setPlacement(envPV);
return det;
}
// helper function to build module with or w/o wrapper
std::tuple<Volume, int, double, double> build_shashlik(Detector &desc, xml::Collection_t &plm, SensitiveDetector &sens)
{
auto mod = plm.child(_Unicode(module));
// a modular volume
std::string shape = dd4hep::getAttrOrDefault(mod, _Unicode(shape), "square");
std::transform(shape.begin(), shape.end(), shape.begin(), [](char c) { return std::tolower(c); });
int nsides = 4;
if (shape == "hexagon") {
nsides = 6;
} else if (shape != "square") {
std::cerr << "ShashlikCalorimeter Error: Unsupported shape of module " << shape
<< ". Please choose from (square, hexagon). Proceed with square shape." << std::endl;
}
double slen = mod.attr<double>(_Unicode(side_length));
double rmax = slen/2./std::sin(M_PI/nsides);
Layering layering(mod);
auto len = layering.totalThickness();
// wrapper info
PolyhedraRegular mpoly(nsides, 0., rmax, len);
Volume mvol("shashlik_module_vol", mpoly, desc.air());
mvol.setVisAttributes(desc.visAttributes(dd4hep::getAttrOrDefault(mod, _Unicode(vis), "GreenVis")));
double gap = 0.;
Volume wvol("shashlik_wrapper_vol");
if (plm.hasChild(_Unicode(wrapper))) {
auto wrap = plm.child(_Unicode(wrapper));
gap = wrap.attr<double>(_Unicode(thickness));
if (gap > 1e-6*mm) {
wvol.setSolid(PolyhedraRegular(nsides, 0., rmax + gap, len));
wvol.setMaterial(desc.material(wrap.attr<std::string>(_Unicode(material))));
wvol.setVisAttributes(desc.visAttributes(dd4hep::getAttrOrDefault(wrap, _Unicode(vis), "WhiteVis")));
wvol.placeVolume(mvol, Position{0., 0., 0.});
}
}
// layer start point
double lz = -len/2.;
int lnum = 1;
// Loop over the sets of layer elements in the detector.
for (xml_coll_t li(mod, _U(layer)); li; ++li) {
int repeat = li.attr<int>(_Unicode(repeat));
// Loop over number of repeats for this layer.
for (int j = 0; j < repeat; j++) {
std::string lname = Form("layer%d", lnum);
double lthick = layering.layer(lnum - 1)->thickness(); // Layer's thickness.
PolyhedraRegular lpoly(nsides, 0., rmax, lthick);
Volume lvol(lname, lpoly, desc.air());
// Loop over the sublayers or slices for this layer.
int snum = 1;
double sz = -lthick/2.;
for (xml_coll_t si(li, _U(slice)); si; ++si) {
std::string sname = Form("slice%d", snum);
double sthick = si.attr<double>(_Unicode(thickness));
PolyhedraRegular spoly(nsides, 0., rmax, sthick);
Volume svol(sname, spoly, desc.material(si.attr<std::string>(_Unicode(material))));
std::string issens = dd4hep::getAttrOrDefault(si, _Unicode(sensitive), "no");
std::transform(issens.begin(), issens.end(), issens.begin(), [](char c) { return std::tolower(c); });
if ((issens == "yes") || (issens == "y") || (issens == "true")) {
svol.setSensitiveDetector(sens);
}
svol.setAttributes(desc,
dd4hep::getAttrOrDefault(si, _Unicode(region), ""),
dd4hep::getAttrOrDefault(si, _Unicode(limits), ""),
dd4hep::getAttrOrDefault(si, _Unicode(vis), "InvisibleNoDaughters"));
// Slice placement.
auto slicePV = lvol.placeVolume(svol, Position(0, 0, sz + sthick/2.));
slicePV.addPhysVolID("slice", snum++);
// Increment Z position of slice.
sz += sthick;
}
// Set region, limitset, and vis of layer.
lvol.setAttributes(desc,
dd4hep::getAttrOrDefault(li, _Unicode(region), ""),
dd4hep::getAttrOrDefault(li, _Unicode(limits), ""),
dd4hep::getAttrOrDefault(li, _Unicode(vis), "InvisibleNoDaughters"));
auto layerPV = mvol.placeVolume(lvol, Position(0, 0, lz + lthick/2));
layerPV.addPhysVolID("layer", lnum++);
// Increment to next layer Z position.
lz += lthick;
}
}
if (gap > 1e-6*mm) {
return std::make_tuple(wvol, nsides, 2.*std::sin(M_PI/nsides)*(rmax + gap), len);
} else {
return std::make_tuple(mvol, nsides, slen, len);
}
}
// place disk of modules
static void add_disk_shashlik(Detector& desc, Assembly &env, xml::Collection_t &plm, SensitiveDetector &sens, int sid)
{
auto [mvol, nsides, sidelen, len] = build_shashlik(desc, plm, sens);
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 = (nsides == 6) ? athena::geo::fillHexagons({0., 0.}, sidelen, rmin, rmax, phimin, phimax)
: athena::geo::fillSquares({0., 0.}, sidelen*1.414, 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;
for (auto &p : points) {
Transform3D tr = RotationZYX(rot.z(), rot.y(), rot.x())
* Translation3D(pos.x() + p.x(), pos.y() + p.y(), pos.z() + len/2.)
* RotationZ((nsides == 4) ? 45*degree : 0);
auto modPV = env.placeVolume(mvol, tr);
modPV.addPhysVolID("sector", sector_id).addPhysVolID("module", id_begin + mid++);
}
}
DECLARE_DETELEMENT(ShashlikCalorimeter, create_detector)