From b4e815b397e3370f23c709f0fdbe1bab40c4ba03 Mon Sep 17 00:00:00 2001
From: Chao Peng <cpeng@anl.gov>
Date: Mon, 12 Jul 2021 04:48:47 +0000
Subject: [PATCH] Interlayer barrel ecal
---
athena.xml | 33 +-
bin/make_dawn_views | 14 +-
compact/definitions.xml | 4 +-
compact/ecal.xml | 15 -
compact/ecal_barrel_interlayers.xml | 158 +++++++++
src/BarrelCalorimeterInterlayers_geo.cpp | 397 +++++++++++++++++++++++
6 files changed, 585 insertions(+), 36 deletions(-)
delete mode 100644 compact/ecal.xml
create mode 100644 compact/ecal_barrel_interlayers.xml
create mode 100644 src/BarrelCalorimeterInterlayers_geo.cpp
diff --git a/athena.xml b/athena.xml
index 25eff83a..d54c46d6 100644
--- a/athena.xml
+++ b/athena.xml
@@ -1,5 +1,5 @@
-<lccdd xmlns:compact="http://www.lcsim.org/schemas/compact/1.0"
- xmlns:xs="http://www.w3.org/2001/XMLSchema"
+<lccdd xmlns:compact="http://www.lcsim.org/schemas/compact/1.0"
+ xmlns:xs="http://www.w3.org/2001/XMLSchema"
xs:noNamespaceSchemaLocation="http://www.lcsim.org/schemas/compact/1.0/compact.xsd">
<!-- Some information about detector -->
@@ -8,7 +8,7 @@
url="https://eicweb.phy.anl.gov/EIC/detectors/athena.git"
status="development"
version="v1 2021-03-16">
- <comment>Athena Detector</comment>
+ <comment>Athena Detector</comment>
</info>
<define>
@@ -112,15 +112,15 @@
<include ref="ip6/beampipe.xml"/>
<detectors>
- <detector id="VertexBarrelSubAssembly_ID"
- name="VertexBarrelSubAssembly"
- type="DD4hep_SubdetectorAssembly"
+ <detector id="VertexBarrelSubAssembly_ID"
+ name="VertexBarrelSubAssembly"
+ type="DD4hep_SubdetectorAssembly"
vis="TrackerSubAssemblyVis">
<composite name="VertexBarrel" />
</detector>
- <detector id="VertexEndcapSubAssembly_ID"
- name="VertexEndcapSubAssembly"
- type="DD4hep_SubdetectorAssembly"
+ <detector id="VertexEndcapSubAssembly_ID"
+ name="VertexEndcapSubAssembly"
+ type="DD4hep_SubdetectorAssembly"
vis="TrackerSubAssemblyVis">
<composite name="VertexEndcapN" />
<composite name="VertexEndcapP" />
@@ -142,7 +142,7 @@
<composite name="TrackerEndcapP_Outer"/>
<composite name="TrackerEndcapN_Outer"/>
</detector>
-
+
<detector id="TOFSubAssembly_ID"
name="TOFSubAssembly"
type="DD4hep_SubdetectorAssembly"
@@ -188,14 +188,21 @@
<include ref="compact/central_tracker.xml"/>
<include ref="compact/tof_barrel.xml"/>
-
+
<!--include ref="compact/rwell_tracker_barrel.xml"/-->
<include ref="compact/cb_DIRC.xml"/>
<!-- When changing magnet, also select dimensions in definitions.xml. -->
<include ref="compact/solenoid.xml"/>
- <include ref="compact/ecal.xml"/>
+ <include ref="compact/ci_ecal.xml"/>
+ <!--<include ref="compact/ci_ecal_shashlik.xml"/>-->
+ <!--<include ref="compact/ce_ecal.xml"/>-->
+ <include ref="compact/ce_ecal_crystal_glass.xml"/>
+ <!-- <include ref="compact/ecal_barrel.xml"/> -->
+ <!-- <include ref="compact/ecal_barrel_hybrid.xml"/> -->
+ <include ref="compact/ecal_barrel_interlayers.xml"/>
+
<include ref="compact/hcal.xml"/>
<!--include ref="compact/ce_GEM.xml"/-->
<!--include ref="compact/gem_tracker_endcap.xml"/-->
@@ -213,7 +220,7 @@
<include ref="compact/mm_tracker_barrel.xml"/>
<include ref="compact/cb_VTX_Barrel.xml"/>
-->
-
+
<readouts>
</readouts>
diff --git a/bin/make_dawn_views b/bin/make_dawn_views
index fe13723b..ea438d27 100755
--- a/bin/make_dawn_views
+++ b/bin/make_dawn_views
@@ -104,13 +104,13 @@ if [ "${DETECTOR_ONLY}" -eq "1" ] ; then
-o derp.root -n 1 \
--ui csh --vis -b -m macro/dawn_picture.mac &
-sleep 10
-echo "sleeping 20 secs .. "
-sleep 10
-echo "sleeping 10 secs "
-sleep 5
-echo "sleeping 5 secs "
-sleep 5
+timeout=200
+while [ $timeout -ge 0 ] && [ ! -f "g4_0000.prim" ]
+do
+ echo "terminating in $timeout secs ..."
+ sleep 5
+ ((timeout=timeout-5))
+done
kill %1
else
diff --git a/compact/definitions.xml b/compact/definitions.xml
index ac42e8b4..27d5db26 100644
--- a/compact/definitions.xml
+++ b/compact/definitions.xml
@@ -180,8 +180,9 @@
EndcapP subassembly ID: 102
EndcapN subassembly ID: 103
Crystal subassembly ID: 104
+ Barrel2 subassembly ID: 105
- Unused IDs: 105-109
+ Unused IDs: 106-109
</comment>
<constant name="ECalSubAssembly_ID" value="100"/>
@@ -189,6 +190,7 @@
<constant name="ECalEndcapP_ID" value="102"/>
<constant name="ECalEndcapN_ID" value="103"/>
<constant name="CrystalEndcap_ID" value="104"/>
+ <constant name="ECalBarrel2_ID" value="105"/>
<comment>
=====================================
diff --git a/compact/ecal.xml b/compact/ecal.xml
deleted file mode 100644
index 60704e6c..00000000
--- a/compact/ecal.xml
+++ /dev/null
@@ -1,15 +0,0 @@
-<lccdd>
-
- <comment>barrel ecal</comment>
- <!-- <include ref="ecal_barrel.xml"/> -->
- <include ref="ecal_barrel_hybrid.xml"/>
-
- <comment>electron endcap ecal</comment>
- <!--<include ref="ce_ecal.xml"/>-->
- <include ref="ce_ecal_crystal_glass.xml"/>
-
- <comment>hadron endcap ecal</comment>
- <include ref="ci_ecal.xml"/>
- <!--<include ref="ci_ecal_shashlik.xml"/>-->
-
-</lccdd>
diff --git a/compact/ecal_barrel_interlayers.xml b/compact/ecal_barrel_interlayers.xml
new file mode 100644
index 00000000..8d1c6910
--- /dev/null
+++ b/compact/ecal_barrel_interlayers.xml
@@ -0,0 +1,158 @@
+<lccdd>
+
+ <display>
+ <vis name="EcalBarrelEnvelope_vis" alpha="0.9" r="0.99" g="0.5" b="0" showDaughters="true" visible="false" />
+ <vis name="EcalBarrelStave_vis" alpha="0.9" r="0.99" g="0.5" b="0" showDaughters="true" visible="false" />
+ </display>
+ <define>
+ <comment>
+ ---------------------------------------
+ EM Calorimeter Parameters with AstroPix
+ ---------------------------------------
+ </comment>
+ <constant name="EcalBarrel_Support_thickness" value="5*cm"/>
+ <constant name="EcalBarrel_SiliconThickness" value="500*um"/>
+ <constant name="EcalBarrel_ElectronicsThickness" value="150*um"/>
+ <constant name="EcalBarrel_CopperThickness" value="100*um"/>
+ <constant name="EcalBarrel_KaptonThickness" value="200*um"/>
+ <constant name="EcalBarrel_EpoxyThickness" value="100*um"/>
+ <constant name="EcalBarrel_CarbonThickness" value="0.5*mm"/>
+ <constant name="EcalBarrel_CarbonSpacerWidth" value="4*mm"/>
+ <constant name="EcalBarrel_LayerSpacing" value="10.0*mm"/>
+ <constant name="EcalBarrel_FiberRadius" value="0.5*mm"/>
+ <constant name="EcalBarrel_FiberXSpacing" value="1.34*mm"/>
+ <constant name="EcalBarrel_FiberZSpacing" value="1.22*mm"/>
+ <constant name="EcalBarrel_SpaceBetween" value="1*mm"/>
+ <comment>
+ For Pb/SiFi (GlueX): X0 ~ 1.45 cm
+ For W/SiFi (sPHENIX): X0 ~ 0.7 cm (but different fiber orientation)
+ </comment>
+ <constant name="EcalBarrel_RadiatorThickness" value="EcalBarrel_FiberZSpacing*12"/>
+ <constant name="EcalBarrel_ModRepeat" value="CaloSides"/>
+ <constant name="EcalBarrel_ModLength" value="0.5*m"/>
+ <constant name="EcalBarrel_ModWidth" value="0.5*m"/>
+ <constant name="EcalBarrel_AvailThickness" value="EcalBarrel_TotalThickness-EcalBarrel_Support_thickness"/>
+ <constant name="EcalBarrel_ImagingLayerThickness"
+ value="EcalBarrel_SiliconThickness
+ + EcalBarrel_ElectronicsThickness
+ + EcalBarrel_CopperThickness
+ + EcalBarrel_KaptonThickness
+ + EcalBarrel_EpoxyThickness
+ + EcalBarrel_CarbonThickness
+ + EcalBarrel_LayerSpacing"/>
+
+ <constant name="EcalBarrelImagingLayers_nMax" value="6"/>
+ <constant name="EcalBarrelImagingLayers_num"
+ value="min(EcalBarrelImagingLayers_nMax,
+ floor(EcalBarrel_AvailThickness/
+ (EcalBarrel_ImagingLayerThickness+EcalBarrel_RadiatorThickness+EcalBarrel_SpaceBetween)))"/>
+ <constant name="EcalBarrel_FiberLayerThickness_max"
+ value="max(0, EcalBarrel_AvailThickness-
+ (EcalBarrelImagingLayers_num*EcalBarrel_ImagingLayerThickness))"/>
+ <constant name="EcalBarrel_FiberLayerThickness"
+ value="min(EcalBarrel_FiberLayerThickness_max, EcalBarrel_FiberZSpacing*12*15)"/>
+ </define>
+
+ <limits>
+ </limits>
+
+ <regions>
+ </regions>
+
+ <display>
+ </display>
+
+ <detectors>
+
+ <comment>
+ ---------------------------------------
+ Imaging Layers of Barrel EM Calorimeter
+ ---------------------------------------
+ Silicon (Astropix) readout layers for imaging 3d showers
+ </comment>
+ <detector
+ id="ECalBarrel_ID"
+ name="EcalBarrelImaging"
+ type="athena_EcalBarrelInterlayers"
+ readout="EcalBarrelHits"
+ calorimeterType="EM_BARREL"
+ vis="EcalBarrelEnvelope_vis"
+ offset="EcalBarrel_offset">
+ <dimensions numsides="EcalBarrel_ModRepeat"
+ rmin="EcalBarrel_rmin"
+ z="EcalBarrel_length"/>
+ <staves vis="EcalBarrelStave_vis"/>
+ <layer repeat="EcalBarrelImagingLayers_num" vis="AnlBlue"
+ space_between="EcalBarrel_RadiatorThickness + EcalBarrel_SpaceBetween"
+ space_before="0.*cm">
+ <slice material="Silicon" thickness="EcalBarrel_SiliconThickness" sensitive="yes" limits="cal_limits" vis="AnlGray"/>
+ <slice material="Silicon" thickness="EcalBarrel_ElectronicsThickness" vis="AnlGold"/>
+ <slice material="Copper" thickness="EcalBarrel_CopperThickness" vis="AnlGray"/>
+ <slice material="Kapton" thickness="EcalBarrel_KaptonThickness" vis="AnlGold"/>
+ <slice material="Epoxy" thickness="EcalBarrel_EpoxyThickness" vis="AnlGray"/>
+ <slice material="CarbonFiber" thickness="EcalBarrel_CarbonThickness" vis="AnlGold"/>
+ <slice material="Air" thickness="EcalBarrel_LayerSpacing " vis="AnlGold"/>
+ </layer>
+ </detector>
+
+ <comment>
+ ---------------------------------------
+ Pb/ScFi Layers of Barrel EM Calorimeter
+ ---------------------------------------
+ </comment>
+ <detector
+ id="ECalBarrel2_ID"
+ name="EcalBarrelScFi"
+ type="athena_EcalBarrelInterlayers"
+ readout="EcalBarrelScFiHits"
+ calorimeterType="EM_BARREL"
+ vis="EcalBarrelEnvelope_vis"
+ offset="EcalBarrel_offset">
+ <dimensions numsides="EcalBarrel_ModRepeat"
+ rmin="EcalBarrel_rmin"
+ z="EcalBarrel_length"/>
+ <staves vis="EcalBarrelStave_vis">
+ <support inside="true" material="Steel235" vis="AnlOrange"
+ thickness="EcalBarrel_Support_thickness"
+ n_beams="3" grid_size="25.0*cm" >
+ </support>
+ </staves>
+ <layer repeat="EcalBarrelImagingLayers_num-1" vis="AnlBlue"
+ space_between="EcalBarrel_ImagingLayerThickness + EcalBarrel_SpaceBetween"
+ space_before="EcalBarrel_ImagingLayerThickness + EcalBarrel_SpaceBetween/2.">
+ <slice material="Lead" thickness="EcalBarrel_RadiatorThickness" vis="EcalBarrelFibersVis">
+ <fiber material="PlasticScint"
+ sensitive="yes"
+ vis="EcalBarrelFiberVis"
+ radius="EcalBarrel_FiberRadius"
+ spacing_x="EcalBarrel_FiberXSpacing"
+ spacing_z="EcalBarrel_FiberZSpacing"/>
+ </slice>
+ </layer>
+ <layer repeat="1" vis="AnlBlue"
+ space_before="EcalBarrel_ImagingLayerThickness + EcalBarrel_SpaceBetween">
+ <slice material="Lead" thickness="EcalBarrel_FiberLayerThickness" vis="EcalBarrelFiberLayerVis">
+ <fiber material="PlasticScint"
+ sensitive="yes"
+ vis="EcalBarrelFiberVis"
+ radius="EcalBarrel_FiberRadius"
+ spacing_x="EcalBarrel_FiberXSpacing"
+ spacing_z="EcalBarrel_FiberZSpacing">
+ </fiber>
+ </slice>
+ </layer>
+ </detector>
+ </detectors>
+
+ <readouts>
+ <readout name="EcalBarrelHits">
+ <segmentation type="CartesianGridXY" grid_size_x="0.5 * mm" grid_size_y="0.5 * mm"/>
+ <id>system:8,module:8,layer:8,slice:8,x:32:-16,y:-16</id>
+ </readout>
+ <readout name="EcalBarrelScFiHits">
+ <segmentation type="NoSegmentation"/>
+ <id>system:8,module:8,layer:8,slice:8,grid:16,fiber:16</id>
+ </readout>
+ </readouts>
+
+</lccdd>
diff --git a/src/BarrelCalorimeterInterlayers_geo.cpp b/src/BarrelCalorimeterInterlayers_geo.cpp
new file mode 100644
index 00000000..ec888447
--- /dev/null
+++ b/src/BarrelCalorimeterInterlayers_geo.cpp
@@ -0,0 +1,397 @@
+// Detector plugin to support a hybrid central barrel calorimeter
+// The detector consists of interlayers of Pb/ScFi (segmentation in global r, phi) and W/Si (segmentation in local x, y)
+// Assembly is used as the envelope so two different detectors can be interlayered with each other
+//
+//
+// Implementation of the Sci Fiber geometry: M. Żurek 06/19/2021
+// Support interlayers between multiple detectors: C. Peng 07/09/2021
+
+
+#include "DD4hep/DetFactoryHelper.h"
+#include "XML/Layering.h"
+#include "Math/Point2D.h"
+#include "TGeoPolygon.h"
+
+using namespace std;
+using namespace dd4hep;
+using namespace dd4hep::detail;
+
+typedef ROOT::Math::XYPoint Point;
+// headers for helper functions
+vector<Point> _fiberPositions(double radius, double x_spacing, double z_spacing,
+ double x, double z, double phi, double spacing_tol = 1e-2);
+std::pair<int, int> _getNdivisions(double x, double z, double dx, double dz);
+vector<tuple<int, Point, Point, Point, Point>> _gridPoints(int div_x, int div_z, double x, double z, double phi);
+
+
+// barrel ecal layers contained in an assembly
+static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector sens) {
+ Layering layering (e);
+ xml_det_t x_det = e;
+ Material air = description.air();
+ int det_id = x_det.id();
+ string det_name = x_det.nameStr();
+ xml_comp_t x_staves = x_det.staves();
+ double offset = x_det.attr<double>(_Unicode(offset));
+ xml_comp_t x_dim = x_det.dimensions();
+ int nsides = x_dim.numsides();
+ double inner_r = x_dim.rmin();
+ double dphi = (2*M_PI/nsides);
+ double hphi = dphi/2;
+
+ DetElement sdet (det_name, det_id);
+ Volume motherVol = description.pickMotherVolume(sdet);
+
+ Assembly envelope (det_name);
+ Transform3D tr = Translation3D(0, 0, offset) * RotationZ(hphi);
+ PlacedVolume env_phv = motherVol.placeVolume(envelope, tr);
+ sens.setType("calorimeter");
+
+ env_phv.addPhysVolID("system",det_id);
+ sdet.setPlacement(env_phv);
+
+ // build a single stave
+ DetElement stave_det("stave0", det_id);
+ Assembly mod_vol("stave");
+
+ // keep tracking of the total thickness
+ double l_pos_z = inner_r;
+ { // ===== buildBarrelStave(description, sens, module_volume) =====
+ // Parameters for computing the layer X dimension:
+ double tan_hphi = std::tan(hphi);
+ double l_dim_y = x_dim.z()/2.;
+
+ // Loop over the sets of layer elements in the detector.
+ int l_num = 1;
+ for(xml_coll_t li(x_det, _U(layer)); li; ++li) {
+ xml_comp_t x_layer = li;
+ int repeat = x_layer.repeat();
+ double l_space_between = dd4hep::getAttrOrDefault(x_layer, _Unicode(space_between), 0.);
+ double l_space_before = dd4hep::getAttrOrDefault(x_layer, _Unicode(space_before), 0.);
+ l_pos_z += l_space_before;
+ // Loop over number of repeats for this layer.
+ for (int j = 0; j < repeat; j++) {
+ string l_name = Form("layer%d", l_num);
+ double l_thickness = layering.layer(l_num - 1)->thickness(); // Layer's thickness.
+ double l_dim_x = tan_hphi* l_pos_z;
+ l_pos_z += l_thickness;
+
+ Position l_pos(0, 0, l_pos_z - l_thickness/2.); // Position of the layer.
+ double l_trd_x1 = l_dim_x;
+ double l_trd_x2 = l_dim_x + l_thickness*tan_hphi;
+ double l_trd_y1 = l_dim_y;
+ double l_trd_y2 = l_trd_y1;
+ double l_trd_z = l_thickness/2;
+ Trapezoid l_shape(l_trd_x1, l_trd_x2, l_trd_y1, l_trd_y2, l_trd_z);
+ Volume l_vol(l_name, l_shape, air);
+ DetElement layer(stave_det, l_name, det_id);
+
+ // Loop over the sublayers or slices for this layer.
+ int s_num = 1;
+ double s_pos_z = -(l_thickness / 2.);
+ for(xml_coll_t si(x_layer,_U(slice)); si; ++si) {
+ xml_comp_t x_slice = si;
+ string s_name = Form("slice%d", s_num);
+ double s_thick = x_slice.thickness();
+ double s_trd_x1 = l_dim_x + (s_pos_z + l_thickness/2)*tan_hphi;
+ double s_trd_x2 = l_dim_x + (s_pos_z + l_thickness/2 + s_thick)*tan_hphi;
+ double s_trd_y1 = l_trd_y1;
+ double s_trd_y2 = s_trd_y1;
+ double s_trd_z = s_thick/2.;
+ Trapezoid s_shape(s_trd_x1, s_trd_x2, s_trd_y1, s_trd_y2, s_trd_z);
+ Volume s_vol(s_name, s_shape, description.material(x_slice.materialStr()));
+ DetElement slice(layer, s_name, det_id);
+
+ // build fibers
+ if (x_slice.hasChild("fiber")) {
+ xml_comp_t x_fiber = x_slice.child(_Unicode(fiber));
+ double f_radius = getAttrOrDefault(x_fiber, _U(radius), 0.1 * cm);
+ double f_spacing_x = getAttrOrDefault(x_fiber, _Unicode(spacing_x), 0.122 * cm);
+ double f_spacing_z = getAttrOrDefault(x_fiber, _Unicode(spacing_z), 0.134 * cm);
+ std::string f_id_grid = getAttrOrDefault(x_fiber, _Unicode(identifier_grid), "grid");
+ std::string f_id_fiber = getAttrOrDefault(x_fiber, _Unicode(identifier_fiber), "fiber");
+
+ // Calculate fiber positions inside the slice
+ vector<Point> f_pos = _fiberPositions(f_radius, f_spacing_x, f_spacing_z, s_trd_x1, s_thick, hphi);
+ // Sort fiber IDs fo better organization
+ sort(f_pos.begin(), f_pos.end(),
+ [](const Point &p1, const Point &p2) {
+ if (p1.y() == p2.y()) { return p1.x() < p2.x(); }
+ return p1.y() < p2.y();
+ });
+
+ Tube f_tube(0, f_radius, l_dim_y);
+
+ // Set up the readout grid for the fiber layers
+ // Trapezoid is divided into segments with equal dz and equal number of divisions in x
+ // Every segment is a polygon that can be attached later to the lightguide
+ // The grid size is assumed to be ~2x2 cm (starting values). This is to be larger than
+ // SiPM chip (for GlueX 13mmx13mm: 4x4 grid 3mmx3mm with 3600 50×50 μm pixels each)
+ // See, e.g., https://arxiv.org/abs/1801.03088 Fig. 2d
+
+ // Calculate number of divisions
+ auto grid_div = _getNdivisions(s_trd_x1, s_thick, 2.0*cm, 2.0*cm);
+ // Calculate polygonal grid coordinates (vertices)
+ auto grid_vtx = _gridPoints(grid_div.first, grid_div.second, s_trd_x1, s_thick, hphi);
+
+ vector<int> f_id_count(grid_div.first*grid_div.second, 0);
+ for (auto &p : f_pos) {
+ int f_grid_id = -1;
+ int f_id = -1;
+ // Check to which grid fiber belongs to
+ for (auto &poly_vtx : grid_vtx) {
+ auto [grid_id, vtx_a, vtx_b, vtx_c, vtx_d] = poly_vtx;
+ double poly_x[4] = {vtx_a.x(), vtx_b.x(), vtx_c.x(), vtx_d.x()};
+ double poly_y[4] = {vtx_a.y(), vtx_b.y(), vtx_c.y(), vtx_d.y()};
+ double f_xy[2] = {p.x(), p.y()};
+
+ TGeoPolygon poly(4);
+ poly.SetXY(poly_x,poly_y);
+ poly.FinishPolygon();
+
+ if(poly.Contains(f_xy)) {
+ f_grid_id = grid_id;
+ f_id = f_id_count[grid_id];
+ f_id_count[grid_id]++;
+ }
+ }
+
+ string f_name = "fiber" + to_string(f_grid_id) + "_" + to_string(f_id);
+ Volume f_vol(f_name, f_tube, description.material(x_fiber.materialStr()));
+ DetElement fiber(slice, f_name, det_id);
+ if ( x_fiber.isSensitive() ) {
+ f_vol.setSensitiveDetector(sens);
+ }
+ fiber.setAttributes(description,f_vol,x_fiber.regionStr(),x_fiber.limitsStr(),x_fiber.visStr());
+
+ // Fiber placement
+ Transform3D f_tr(RotationZYX(0,0,M_PI*0.5),Position(p.x(), 0 ,p.y()));
+ PlacedVolume fiber_phv = s_vol.placeVolume(f_vol, f_tr);
+ fiber_phv.addPhysVolID(f_id_grid, f_grid_id + 1).addPhysVolID(f_id_fiber, f_id + 1);
+ fiber.setPlacement(fiber_phv);
+ }
+ }
+
+
+ if ( x_slice.isSensitive() ) {
+ s_vol.setSensitiveDetector(sens);
+ }
+ s_vol.setAttributes(description, x_slice.regionStr(), x_slice.limitsStr(), x_slice.visStr());
+
+ // Slice placement.
+ PlacedVolume slice_phv = l_vol.placeVolume(s_vol, Position(0, 0, s_pos_z + s_thick/2));
+ slice_phv.addPhysVolID("slice", s_num);
+ slice.setPlacement(slice_phv);
+ // Increment Z position of slice.
+ s_pos_z += s_thick;
+ ++s_num;
+ }
+
+ // Set region, limitset, and vis of layer.
+ l_vol.setAttributes(description, x_layer.regionStr(), x_layer.limitsStr(), x_layer.visStr());
+
+ PlacedVolume layer_phv = mod_vol.placeVolume(l_vol, l_pos);
+ layer_phv.addPhysVolID("layer", l_num);
+ layer.setPlacement(layer_phv);
+ // Increment to next layer Z position. Do not add space_between for the last layer
+ if (j < repeat - 1) {
+ l_pos_z += l_space_between;
+ }
+ ++l_num;
+ }
+ }
+ }
+ // Phi start for a stave.
+ double phi = M_PI / nsides;
+ // Create nsides staves.
+ for (int i = 0; i < nsides; i++, phi -= dphi) { // i is module number
+ // Compute the stave position
+ Transform3D tr(RotationZYX(0, phi, M_PI*0.5), Translation3D(0, 0, 0));
+ PlacedVolume pv = envelope.placeVolume(mod_vol, tr);
+ pv.addPhysVolID("module", i + 1);
+ DetElement sd = (i == 0) ? stave_det : stave_det.clone(Form("stave%d", i));
+ sd.setPlacement(pv);
+ sdet.add(sd);
+ }
+
+ Solid support_frame_s;
+ // optional stave support
+ if (x_staves.hasChild("support")) {
+ xml_comp_t x_support = x_staves.child(_U(support));
+ double support_thickness = getAttrOrDefault(x_support, _U(thickness), 5.0 * cm);
+ double trd_x1_support = (2 * std::tan(hphi) * l_pos_z + support_thickness)/2;
+ // is the support on the inside surface?
+ bool is_inside_support = getAttrOrDefault<bool>(x_support, _Unicode(inside), true);
+ double trd_x1 = std::tan(hphi) * inner_r;
+ double trd_x2 = std::tan(hphi) * (l_pos_z + support_thickness);
+ double trd_y1 = x_dim.z()/2.;
+
+ // 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) * (l_pos_z + beam_thickness)) / 2.;
+ 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.;
+ 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 , trd_y1, support_thickness / 2.0 );
+ Box beam_hori_s(beam_width / 2.0, trd_y1, beam_thickness / 2.0);
+ 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, beam_thickness / 2.0,
+ support_thickness / 2.0 - 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));
+ }
+
+ //l_pos_z += support_thickness;
+ // Set envelope volume attributes.
+ envelope.setAttributes(description, x_det.regionStr(), x_det.limitsStr(), x_det.visStr());
+ return sdet;
+}
+
+
+DECLARE_DETELEMENT(athena_EcalBarrelInterlayers, create_detector)
+// DECLARE_DETELEMENT(athena_EcalBarrelInterlayers, create_detector)
+
+
+// -----------------------------------------------------------------------------
+// helper functions
+// -----------------------------------------------------------------------------
+// Fill fiber lattice into trapezoid starting from position (0,0) in x-z coordinate system
+vector<Point> _fiberPositions(double radius, double x_spacing, double z_spacing,
+ double x, double z, double phi, double spacing_tol) {
+ // z_spacing - distance between fiber layers in z
+ // x_spacing - distance between fiber centers in x
+ // x - half-length of the shorter (bottom) base of the trapezoid
+ // z - height of the trapezoid
+ // phi - angle between z and trapezoid arm
+
+ vector<Point> positions;
+ int z_layers = floor((z/2-radius-spacing_tol)/z_spacing); // number of layers that fit in z/2
+
+ double z_pos = 0.;
+ double x_pos = 0.;
+
+ for(int l = -z_layers; l < z_layers+1; l++) {
+
+ z_pos = l*z_spacing;
+ double x_max = x + (z/2. + z_pos)*tan(phi) - spacing_tol; // calculate max x at particular z_pos
+ (l % 2 == 0) ? x_pos = 0. : x_pos = x_spacing/2; // account for spacing/2 shift
+
+ while(x_pos < (x_max - radius)) {
+ positions.push_back(Point(x_pos,z_pos));
+ if(x_pos != 0.) positions.push_back(Point(-x_pos,z_pos)); // using symmetry around x=0
+ x_pos += x_spacing;
+ }
+ }
+
+ return positions;
+}
+
+// Calculate number of divisions for the readout grid for the fiber layers
+std::pair<int, int> _getNdivisions(double x, double z, double dx, double dz){
+ // x and z defined as in vector<Point> fiberPositions
+ // dx, dz - size of the grid in x and z we want to get close to with the polygons
+ // See also descripltion when the function is called
+
+ double SiPMsize = 13.0*mm;
+ double grid_min = SiPMsize + 3.0*mm;
+
+ if(dz < grid_min) {
+ dz = grid_min;
+ }
+
+ if(dx < grid_min) {
+ dx = grid_min;
+ }
+
+ int nfit_cells_z = floor(z/dz);
+ int n_cells_z = nfit_cells_z;
+
+ if(nfit_cells_z == 0) n_cells_z++;
+
+ int nfit_cells_x = floor((2*x)/dx);
+ int n_cells_x = nfit_cells_x;
+
+ if(nfit_cells_x == 0) n_cells_x++;
+
+ return std::make_pair(n_cells_x, n_cells_z);
+
+}
+
+// Calculate dimensions of the polygonal grid in the cartesian coordinate system x-z
+vector< tuple<int, Point, Point, Point, Point> > _gridPoints(int div_x, int div_z, double x, double z, double phi) {
+ // x, z and phi defined as in vector<Point> fiberPositions
+ // div_x, div_z - number of divisions in x and z
+ double dz = z/div_z;
+
+ std::vector<std::tuple<int, Point, Point, Point, Point>> points;
+
+ for(int iz = 0; iz < div_z + 1; iz++){
+ for(int ix = 0; ix < div_x + 1; ix++){
+ double A_z = -z/2 + iz*dz;
+ double B_z = -z/2 + (iz+1)*dz;
+
+ double len_x_for_z = 2*(x+iz*dz*tan(phi));
+ double len_x_for_z_plus_1 = 2*(x + (iz+1)*dz*tan(phi));
+
+ double dx_for_z = len_x_for_z/div_x;
+ double dx_for_z_plus_1 = len_x_for_z_plus_1/div_x;
+
+ double A_x = -len_x_for_z/2. + ix*dx_for_z;
+ double B_x = -len_x_for_z_plus_1/2. + ix*dx_for_z_plus_1;
+
+ double C_z = B_z;
+ double D_z = A_z;
+ double C_x = B_x + dx_for_z_plus_1;
+ double D_x = A_x + dx_for_z;
+
+ int id = ix + div_x * iz;
+
+ auto A = Point(A_x, A_z);
+ auto B = Point(B_x, B_z);
+ auto C = Point(C_x, C_z);
+ auto D = Point(D_x, D_z);
+
+ // vertex points filled in the clock-wise direction
+ points.push_back(make_tuple(id, A, B, C, D));
+
+ }
+ }
+
+ return points;
+
+}
+
--
GitLab