From 15115f460fb7f1428816369a7bb04e326d7064cf Mon Sep 17 00:00:00 2001
From: Chao Peng <cpeng@anl.gov>
Date: Thu, 1 Jul 2021 02:24:59 +0000
Subject: [PATCH] implement Shashlik Calorimeter
---
athena.xml | 6 -
compact/ci_ecal.xml | 76 ++++++++++++
compact/ci_ecal_shashlik.xml | 76 ++++++++++++
compact/ecal.xml | 84 +------------
src/GaseousRICH_geo.cpp | 2 +-
src/GeometryHelpers.cpp | 134 ++++++++++++++------
src/GeometryHelpers.h | 7 +-
src/HomogeneousCalorimeter_geo.cpp | 28 ++---
src/ShashlikCalorimeter.cpp | 189 +++++++++++++++++++++++++++++
src/ce_MRICH.cpp | 2 +-
10 files changed, 462 insertions(+), 142 deletions(-)
create mode 100644 compact/ci_ecal.xml
create mode 100644 compact/ci_ecal_shashlik.xml
create mode 100644 src/ShashlikCalorimeter.cpp
diff --git a/athena.xml b/athena.xml
index a4c3f87d..ff84f5a1 100644
--- a/athena.xml
+++ b/athena.xml
@@ -108,8 +108,6 @@
</display>
<comment> Include the IP components first </comment>
- <include ref="ip6/forward_ion_beamline.xml"/>
- <include ref="ip6/beampipe.xml"/>
<detectors>
<detector id="VertexBarrelSubAssembly_ID"
@@ -203,10 +201,6 @@
<include ref="compact/forward_trd.xml"/>
<include ref="compact/gaseous_rich.xml"/>
- <include ref="ip6/B0_tracker.xml"/>
- <include ref="ip6/far_forward_offM_tracker.xml"/>
- <include ref="ip6/far_forward_romanpots.xml"/>
- <include ref="ip6/far_forward_detectors.xml"/>
<!--
<include ref="compact/mm_tracker_barrel.xml"/>
diff --git a/compact/ci_ecal.xml b/compact/ci_ecal.xml
new file mode 100644
index 00000000..5ea301e7
--- /dev/null
+++ b/compact/ci_ecal.xml
@@ -0,0 +1,76 @@
+<lccdd>
+ <define>
+ <constant name="EcalEndcapP_rmax" value="Solenoid_rmax "/>
+ </define>
+
+
+ <limits>
+ </limits>
+
+ <regions>
+ </regions>
+
+ <!-- Common Generic visualization attributes -->
+ <comment>Common Generic visualization attributes</comment>
+ <display>
+ </display>
+
+ <detectors>
+
+ <comment>
+ ------------------------------------------
+ Forward (Positive Z) Endcap EM Calorimeter
+ ------------------------------------------
+ A layered EM calorimeter with tungsten and silicon (or scintillator) strips
+ </comment>
+ <detector id="ECalEndcapP_ID"
+ name="EcalEndcapP"
+ reflect="false"
+ type="refdet_PolyhedraEndcapCalorimeter2"
+ readout="EcalEndcapPHits"
+ vis="EcalEndcapVis"
+ calorimeterType="EM_ENDCAP" >
+ <position x="0" y="0" z="-0"/>
+ <dimensions
+ numsides="CaloSides"
+ zmin="EcalEndcapP_zmin"
+ rmin="EcalEndcapP_rmin"
+ rmax="EcalEndcapP_rmax " />
+ <layer repeat="EcalEndcapPLayer1_NRepeat">
+ <slice material="Silicon" thickness="EcalSiliconThickness" sensitive="yes" limits="cal_limits"/>
+ <slice material="Copper" thickness="EcalCopperThickness"/>
+ <slice material="Kapton" thickness="EcalKaptonThickness"/>
+ <slice material="Air" thickness="EcalAir1Thickness"/>
+ </layer>
+ <layer repeat="EcalEndcapPLayer2_NRepeat">
+ <slice material="TungstenDens24" thickness="EcalThinTungstenThickness"/>
+ <slice material="Air" thickness="EcalAir2Thickness"/>
+ <slice material="Silicon" thickness="EcalSiliconThickness" sensitive="yes" limits="cal_limits"/>
+ <slice material="Copper" thickness="EcalCopperThickness"/>
+ <slice material="Kapton" thickness="EcalKaptonThickness"/>
+ <slice material="Air" thickness="EcalAir1Thickness"/>
+ </layer>
+ <layer repeat="EcalEndcapPLayer3_NRepeat">
+ <slice material="TungstenDens24" thickness="EcalThickTungstenThickness"/>
+ <slice material="Air" thickness="EcalAir2Thickness"/>
+ <slice material="Silicon" thickness="EcalSiliconThickness" sensitive="yes" limits="cal_limits"/>
+ <slice material="Copper" thickness="EcalCopperThickness"/>
+ <slice material="Kapton" thickness="EcalKaptonThickness"/>
+ <slice material="Air" thickness="EcalAir1Thickness"/>
+ </layer>
+ </detector>
+ </detectors>
+
+ <!-- Definition of the readout segmentation/definition -->
+ <readouts>
+ <readout name="EcalEndcapPHits">
+ <segmentation type="CartesianGridXY" grid_size_x="3.5 * mm" grid_size_y="3.5 * mm"/>
+ <id>system:8,barrel:3,module:4,layer:6,slice:5,x:32:-16,y:-16</id>
+ </readout>
+ </readouts>
+
+ <plugins>
+ </plugins>
+
+
+</lccdd>
diff --git a/compact/ci_ecal_shashlik.xml b/compact/ci_ecal_shashlik.xml
new file mode 100644
index 00000000..171c59a4
--- /dev/null
+++ b/compact/ci_ecal_shashlik.xml
@@ -0,0 +1,76 @@
+<lccdd>
+ <define>
+ <constant name="EcalEndcapP_rmax" value="Solenoid_rmax "/>
+ </define>
+
+
+ <limits>
+ </limits>
+
+ <regions>
+ </regions>
+
+ <!-- Common Generic visualization attributes -->
+ <comment>Common Generic visualization attributes</comment>
+ <display>
+ </display>
+
+ <detectors>
+
+ <comment>
+ ------------------------------------------
+ Forward (Positive Z) Endcap EM Calorimeter
+ ------------------------------------------
+ An EM calorimeter with shashlik hexagon modules
+ </comment>
+ <detector id="ECalEndcapP_ID"
+ name="EcalEndcapP"
+ type="ShashlikCalorimeter"
+ readout="EcalEndcapPHits">
+ <position x="0" y="0" z="EcalEndcapP_zmin"/>
+ <placements>
+ <disk rmin="EcalEndcapP_rmin" rmax = "EcalEndcapP_rmax" sector="1">
+ <wrapper thickness="2*mm" material="Epoxy" vis="WhiteVis"/>
+ <module shape="square" side_length="50*mm" vis="EcalEndcapVis">
+ <layer repeat="EcalEndcapPLayer1_NRepeat" vis="EcalEndcapVis">
+ <slice material="Silicon" thickness="EcalSiliconThickness" sensitive="yes" limits="cal_limits"/>
+ <slice material="Copper" thickness="EcalCopperThickness"/>
+ <slice material="Kapton" thickness="EcalKaptonThickness"/>
+ <slice material="Air" thickness="EcalAir1Thickness"/>
+ </layer>
+ <layer repeat="EcalEndcapPLayer2_NRepeat">
+ <slice material="TungstenDens24" thickness="EcalThinTungstenThickness"/>
+ <slice material="Air" thickness="EcalAir2Thickness"/>
+ <slice material="Silicon" thickness="EcalSiliconThickness" sensitive="yes" limits="cal_limits"/>
+ <slice material="Copper" thickness="EcalCopperThickness"/>
+ <slice material="Kapton" thickness="EcalKaptonThickness"/>
+ <slice material="Air" thickness="EcalAir1Thickness"/>
+ </layer>
+ <layer repeat="EcalEndcapPLayer3_NRepeat">
+ <slice material="TungstenDens24" thickness="EcalThickTungstenThickness"/>
+ <slice material="Air" thickness="EcalAir2Thickness"/>
+ <slice material="Silicon" thickness="EcalSiliconThickness" sensitive="yes" limits="cal_limits"/>
+ <slice material="Copper" thickness="EcalCopperThickness"/>
+ <slice material="Kapton" thickness="EcalKaptonThickness"/>
+ <slice material="Air" thickness="EcalAir1Thickness"/>
+ </layer>
+ </module>
+ </disk>
+ </placements>
+ </detector>
+ </detectors>
+
+ <!-- Definition of the readout segmentation/definition -->
+ <readouts>
+ <readout name="EcalEndcapPHits">
+ <segmentation type="NoSegmentation"/>
+ <id>system:8,sector:4,module:24,layer:8,slice:8</id>
+ </readout>
+ </readouts>
+
+ <plugins>
+ </plugins>
+
+
+</lccdd>
+
diff --git a/compact/ecal.xml b/compact/ecal.xml
index 4e19bb16..60704e6c 100644
--- a/compact/ecal.xml
+++ b/compact/ecal.xml
@@ -1,87 +1,15 @@
<lccdd>
- <define>
- <constant name="EcalEndcapP_rmax" value="Solenoid_rmax "/>
- </define>
-
-
- <limits>
- </limits>
-
- <regions>
- </regions>
-
- <!-- Common Generic visualization attributes -->
- <comment>Common Generic visualization attributes</comment>
- <display>
- </display>
-
+ <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"/>
- <detectors>
-
- <comment>
- ------------------------------------------
- Forward (Positive Z) Endcap EM Calorimeter
- ------------------------------------------
- A layered EM calorimeter with tungsten and silicon (or scintillator) strips
- </comment>
- <detector id="ECalEndcapP_ID"
- name="EcalEndcapP"
- reflect="false"
- type="refdet_PolyhedraEndcapCalorimeter2"
- readout="EcalEndcapPHits"
- vis="EcalEndcapVis"
- calorimeterType="EM_ENDCAP" >
- <position x="0" y="0" z="-0"/>
- <dimensions
- numsides="CaloSides"
- zmin="EcalEndcapP_zmin"
- rmin="EcalEndcapP_rmin"
- rmax="EcalEndcapP_rmax " />
- <layer repeat="EcalEndcapPLayer1_NRepeat">
- <slice material="Silicon" thickness="EcalSiliconThickness" sensitive="yes" limits="cal_limits"/>
- <slice material="Copper" thickness="EcalCopperThickness"/>
- <slice material="Kapton" thickness="EcalKaptonThickness"/>
- <slice material="Air" thickness="EcalAir1Thickness"/>
- </layer>
- <layer repeat="EcalEndcapPLayer2_NRepeat">
- <slice material="TungstenDens24" thickness="EcalThinTungstenThickness"/>
- <slice material="Air" thickness="EcalAir2Thickness"/>
- <slice material="Silicon" thickness="EcalSiliconThickness" sensitive="yes" limits="cal_limits"/>
- <slice material="Copper" thickness="EcalCopperThickness"/>
- <slice material="Kapton" thickness="EcalKaptonThickness"/>
- <slice material="Air" thickness="EcalAir1Thickness"/>
- </layer>
- <layer repeat="EcalEndcapPLayer3_NRepeat">
- <slice material="TungstenDens24" thickness="EcalThickTungstenThickness"/>
- <slice material="Air" thickness="EcalAir2Thickness"/>
- <slice material="Silicon" thickness="EcalSiliconThickness" sensitive="yes" limits="cal_limits"/>
- <slice material="Copper" thickness="EcalCopperThickness"/>
- <slice material="Kapton" thickness="EcalKaptonThickness"/>
- <slice material="Air" thickness="EcalAir1Thickness"/>
- </layer>
- </detector>
- </detectors>
-
- <!-- Definition of the readout segmentation/definition -->
- <readouts>
- <!--
- <readout name="PlaneTrackerHits">
- <segmentation type="CartesianGridXY" grid_size_x="20.0*mm" grid_size_y="20.0*mm" />
- <id>system:5,module:4,x:32:-16,y:-16</id>
- </readout>
- -->
- <readout name="EcalEndcapPHits">
- <segmentation type="CartesianGridXY" grid_size_x="3.5 * mm" grid_size_y="3.5 * mm"/>
- <id>system:8,barrel:3,module:4,layer:6,slice:5,x:32:-16,y:-16</id>
- </readout>
- </readouts>
- <plugins>
- </plugins>
+ <comment>hadron endcap ecal</comment>
+ <include ref="ci_ecal.xml"/>
+ <!--<include ref="ci_ecal_shashlik.xml"/>-->
</lccdd>
diff --git a/src/GaseousRICH_geo.cpp b/src/GaseousRICH_geo.cpp
index 38a66b8d..ab26484b 100644
--- a/src/GaseousRICH_geo.cpp
+++ b/src/GaseousRICH_geo.cpp
@@ -192,7 +192,7 @@ void build_sensors(Detector &desc, Volume &env, xml::Component plm, const Positi
double rotx = dd4hep::getAttrOrDefault<double>(sec, _Unicode(rotx), 0.);
// fill sensors to the piece
- auto points = ref::utils::fillRectangles({0., 0.}, sx + gap, sy + gap, rmin - gap, rmax + gap, -phiw/2., phiw/2.);
+ 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
diff --git a/src/GeometryHelpers.cpp b/src/GeometryHelpers.cpp
index c37d30bd..725aa893 100644
--- a/src/GeometryHelpers.cpp
+++ b/src/GeometryHelpers.cpp
@@ -1,12 +1,23 @@
#include "GeometryHelpers.h"
// some utility functions that can be shared
-namespace ref::utils {
+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 in_ring(const Point &pt, double sx, double sy, double rmin, double rmax, double phmin, double phmax)
+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;
@@ -27,33 +38,29 @@ inline bool in_ring(const Point &pt, double sx, double sy, double rmin, double r
return true;
}
-// check if a square is overlapped with the others
-inline bool overlap(const Point &pt, double sx, double sy, const std::vector<Point> &pts)
-{
- for (auto &p : pts) {
- auto pn = p - pt;
- if ((std::abs(pn.x()) < (1. - 1e-6)*sx) && (std::abs(pn.y()) < (1. - 1e-6)*sy)) {
- return true;
- }
- }
- return false;
-}
-
// 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)
+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)
{
- // outside of the ring or overlapping
- if (!in_ring(p, sx, sy, rmin, rmax, phmin, phmax) || overlap(p, sx, sy, res)) {
+ // 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;
}
- res.emplace_back(p);
+ bool in_ring = rec_in_ring(p, sx, sy, rmin, rmax, phmin, phmax);
+ if (in_ring) {
+ res.emplace_back(p);
+ }
- // check adjacent squares
- add_rectangle(Point(p.x() + sx, p.y()), res, sx, sy, rmin, rmax, phmin, phmax);
- add_rectangle(Point(p.x() - sx, p.y()), res, sx, sy, rmin, rmax, phmin, phmax);
- add_rectangle(Point(p.x(), p.y() + sy), res, sx, sy, rmin, rmax, phmin, phmax);
- add_rectangle(Point(p.x(), p.y() - sy), res, sx, sy, rmin, rmax, phmin, phmax);
+ // 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
@@ -68,23 +75,76 @@ std::vector<Point> fillRectangles(Point ref, double sx, double sy, double rmin,
// move to center
ref = ref - Point(int(ref.x()/sx)*sx, int(ref.y()/sy)*sy);
- auto find_seed = [] (const Point &ref, int n, double sx, double sy, double rmin, double rmax, double phmin, double phmax) {
- for (int ix = -n; ix < n; ++ix) {
- for (int iy = -n; iy < n; ++iy) {
- Point pt(ref.x() + ix*sx, ref.y() + iy*sy);
- if (in_ring(pt, sx, sy, rmin, rmax, phmin, phmax)) {
- return pt;
- }
- }
+ 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 ref;
- };
+ }
+ 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;
- ref = find_seed(ref, int(rmax/sx) + 2, sx, sy, rmin, rmax, phmin, phmax);
- add_rectangle(ref, res, sx, sy, rmin, rmax, phmin, phmax);
+ add_poly(ref, res, 6, lside, rmin, rmax, phmin, phmax, std::pow(int(rmax/lside) + 1, 2)*2);
return res;
}
-
-} // ref::utils
+} // athena::geo
diff --git a/src/GeometryHelpers.h b/src/GeometryHelpers.h
index 028c474c..5ec55025 100644
--- a/src/GeometryHelpers.h
+++ b/src/GeometryHelpers.h
@@ -3,7 +3,7 @@
#include "Math/Point2D.h"
// some utility functions that can be shared
-namespace ref::utils {
+namespace athena::geo {
typedef ROOT::Math::XYPoint Point;
@@ -17,4 +17,7 @@ inline std::vector<Point> fillSquares(Point ref, double size, double rmin, doubl
return fillRectangles(ref, size, size, rmin, rmax, phmin, phmax);
}
-} // ref::utils
+std::vector<Point> fillHexagons(Point ref, double lside, double rmin, double rmax,
+ double phmin = -M_PI, double phmax = M_PI);
+
+} // athena::geo
diff --git a/src/HomogeneousCalorimeter_geo.cpp b/src/HomogeneousCalorimeter_geo.cpp
index a365801d..ec0b3a62 100644
--- a/src/HomogeneousCalorimeter_geo.cpp
+++ b/src/HomogeneousCalorimeter_geo.cpp
@@ -16,10 +16,10 @@
#include <XML/Helper.h>
#include <iostream>
#include <algorithm>
+#include <tuple>
#include <math.h>
using namespace dd4hep;
-using namespace dd4hep::detail;
/** \addtogroup calorimeters Calorimeters
*/
@@ -170,13 +170,12 @@ static Ref_t create_detector(Detector& desc, xml::Handle_t handle, SensitiveDete
}
// helper function to build module with or w/o wrapper
-Volume build_module(Detector &desc, xml::Collection_t &plm, SensitiveDetector &sens, Position &dim)
+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));
- dim = Position{sx, sy, sz};
Box modShape(sx/2., sy/2., sz/2.);
auto modMat = desc.material(mod.attr<std::string>(_Unicode(material)));
Volume modVol("module_vol", modShape, modMat);
@@ -185,29 +184,27 @@ Volume build_module(Detector &desc, xml::Collection_t &plm, SensitiveDetector &s
// no wrapper
if (!plm.hasChild(_Unicode(wrapper))) {
- return modVol;
+ 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 == 0.) {
- return modVol;
+ 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))));
- dim = Position{sx + thickness, sy + thickness, sz};
- return wrpVol;
+ return std::make_tuple(wrpVol, Position{sx + thickness, sy + thickness, sz});
}
}
// place modules, id must be provided
static void add_individuals(Detector& desc, Assembly &env, xml::Collection_t &plm, SensitiveDetector &sens, int sid)
{
- Position modSize;
- auto modVol = build_module(desc, plm, sens, modSize);
+ auto [modVol, modSize] = build_module(desc, plm, sens);
int sector_id = dd4hep::getAttrOrDefault<int>(plm, _Unicode(sector), sid);
for (xml::Collection_t pl(plm, _Unicode(placement)); pl; ++pl) {
@@ -228,8 +225,7 @@ static void add_individuals(Detector& desc, Assembly &env, xml::Collection_t &pl
// place array of modules
static void add_array(Detector& desc, Assembly &env, xml::Collection_t &plm, SensitiveDetector &sens, int sid)
{
- Position modSize;
- auto modVol = build_module(desc, plm, sens, modSize);
+ 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));
@@ -266,8 +262,7 @@ static void add_array(Detector& desc, Assembly &env, xml::Collection_t &plm, Sen
// place disk of modules
static void add_disk(Detector& desc, Assembly &env, xml::Collection_t &plm, SensitiveDetector &sens, int sid)
{
- Position modSize;
- auto modVol = build_module(desc, plm, sens, modSize);
+ 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));
@@ -275,7 +270,7 @@ static void add_disk(Detector& desc, Assembly &env, xml::Collection_t &plm, Sens
double phimin = dd4hep::getAttrOrDefault<double>(plm, _Unicode(phimin), 0.);
double phimax = dd4hep::getAttrOrDefault<double>(plm, _Unicode(phimax), 2.*M_PI);
- auto points = ref::utils::fillRectangles({0., 0.}, modSize.x(), modSize.y(), rmin, rmax, phimin, phimax);
+ 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));
@@ -291,8 +286,7 @@ static void add_disk(Detector& desc, Assembly &env, xml::Collection_t &plm, Sens
// place lines of modules (anchor point is the 0th module of this line)
static void add_lines(Detector& desc, Assembly &env, xml::Collection_t &plm, SensitiveDetector &sens, int sid)
{
- Position modSize;
- auto modVol = build_module(desc, plm, sens, modSize);
+ 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);
diff --git a/src/ShashlikCalorimeter.cpp b/src/ShashlikCalorimeter.cpp
new file mode 100644
index 00000000..5ce04be8
--- /dev/null
+++ b/src/ShashlikCalorimeter.cpp
@@ -0,0 +1,189 @@
+//==========================================================================
+// 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)
+
diff --git a/src/ce_MRICH.cpp b/src/ce_MRICH.cpp
index ccc2a2b3..5920ddc3 100644
--- a/src/ce_MRICH.cpp
+++ b/src/ce_MRICH.cpp
@@ -82,7 +82,7 @@ void addModules(Volume &mother, xml::DetElement &detElem, Detector &desc, Sensit
modVol.setSensitiveDetector(sens);
// place modules in the sectors (disk)
- auto points = ref::utils::fillSquares({0., 0.}, mWidth + mGap, rmin - mGap, rmax + mGap);
+ auto points = athena::geo::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.;
--
GitLab