diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 9896f2f564d6a31035e4b61c05cb8de6aa3be421..5f0625d6402a1866a20c29535558a48235c51c05 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -192,6 +192,8 @@ overlap_check_geant4:full_fast:
script:
## disable fibers in ECAL for normal overlap check
- sed -i '/<fiber/,+6d' ${DETECTOR_PATH}/compact/ecal_barrel_interlayers.xml
+ ## reduce the number of fibers in Hadron EMCal for overlap check
+ - sed -i 's/radius="EcalEndcapP_FiberRadius"/radius="EcalEndcapP_FiberRadius*10"/' ${DETECTOR_PATH}/compact/ci_ecal_scfi.xml
- python scripts/checkOverlaps.py -c ${DETECTOR_PATH}/athena.xml | tee doc/overlap_check_geant4.out
- echo "$(cat doc/overlap_check_geant4.out | grep GeomVol1002 | wc -l) overlaps..."
- if [[ "$(cat doc/overlap_check_geant4.out | grep GeomVol1002 | wc -l)" -gt "0" ]] ; then echo "Overlaps exist!" && false ; fi
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 5d97ae91e55f6d0654a093c8ee3d5e0faae5f441..5ec35297e96c143e440b4fcb3d9ed3199b0a0e71 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -36,6 +36,7 @@ dd4hep_add_plugin(${a_lib_name} SOURCES
src/HybridCalorimeter_geo.cpp
src/MRich_geo.cpp
src/PolyhedraEndcapCalorimeter2_geo.cpp
+ src/ScFiCalorimeter_geo.cpp
src/ShashlikCalorimeter_geo.cpp
src/SimpleDiskDetector_geo.cpp
src/SolenoidCoil_geo.cpp
diff --git a/bin/generate_prim_file b/bin/generate_prim_file
index 2274f06101c6293258debbdcf18aa4d12deb5116..a0593210bd2e13fd48f68ba165f57bee0e4c18b3 100755
--- a/bin/generate_prim_file
+++ b/bin/generate_prim_file
@@ -54,7 +54,7 @@ parser.add_argument('-t', '--tag', type=str,dest='file_tag',
parser.add_argument('--timeout', type=int,
default=60,
help='Timeout in seconds')
-
+
parser.add_argument('passthrough', nargs='*')
args = parser.parse_args()
@@ -71,6 +71,11 @@ args.out_dir = os.path.abspath(args.out_dir)
args.compact = os.path.abspath(args.compact)
macro = os.path.abspath(macro)
+# adjust fiber radius to reduce the number of fibers
+compact_dir = os.path.dirname(os.path.realpath(args.compact))
+ci_ecal = os.path.join(compact_dir, 'compact', 'ci_ecal_scfi.xml')
+os.system('sed -i \'s/radius=\"EcalEndcapP_FiberRadius\"/radius=\"EcalEndcapP_FiberRadius*10\"/\' {}'.format(ci_ecal))
+
prim_file = 'g4_0000.prim'
dawn_env = os.environ.copy()
dawn_env['DAWN_BATCH'] = 'a'
@@ -138,6 +143,9 @@ for proc in psutil.process_iter():
print('kill {}, generated from {}'.format(pinfo, p.pid))
os.kill(pinfo['pid'], signal.SIGTERM)
+# revert the change
+os.system('sed -i \'s/radius=\"EcalEndcapP_FiberRadius*10\"/radius=\"EcalEndcapP_FiberRadius\"/\' {}'.format(ci_ecal))
+
line = b'stderr outputs:\n'
while line:
print(line.decode('utf-8'), end='')
diff --git a/compact/ce_ecal_crystal_glass.xml b/compact/ce_ecal_crystal_glass.xml
index 2634afa199e7ba06251f6311af02329577e82d7a..cea47efc7029ef807bc8f20ba56e4304645104cf 100644
--- a/compact/ce_ecal_crystal_glass.xml
+++ b/compact/ce_ecal_crystal_glass.xml
@@ -60,32 +60,29 @@
<lines sector="2" mirrorx="true" mirrory="true">
<module sizex="GlassModule_sx" sizey="GlassModule_sy" sizez="GlassModule_sz" material="PbGlass" vis="EcalEndcapNModuleVis"/>
<wrapper thickness="GlassModule_wrap" material="Epoxy" vis="InvisibleWithDaughters"/>
- <line axis="x" x="GlassModule_dx/2." y="GlassModule_dy*1/2." z="GlassModule_z0" begin="13" nmods="13"/>
- <line axis="x" x="GlassModule_dx/2." y="GlassModule_dy*3/2." z="GlassModule_z0" begin="13" nmods="13"/>
- <line axis="x" x="GlassModule_dx/2." y="GlassModule_dy*5/2." z="GlassModule_z0" begin="13" nmods="13"/>
- <line axis="x" x="GlassModule_dx/2." y="GlassModule_dy*7/2." z="GlassModule_z0" begin="12" nmods="14"/>
- <line axis="x" x="GlassModule_dx/2." y="GlassModule_dy*9/2." z="GlassModule_z0" begin="12" nmods="14"/>
- <line axis="x" x="GlassModule_dx/2." y="GlassModule_dy*11/2." z="GlassModule_z0" begin="12" nmods="14"/>
- <line axis="x" x="GlassModule_dx/2." y="GlassModule_dy*13/2." z="GlassModule_z0" begin="11" nmods="15"/>
- <line axis="x" x="GlassModule_dx/2." y="GlassModule_dy*15/2." z="GlassModule_z0" begin="10" nmods="15"/>
- <line axis="x" x="GlassModule_dx/2." y="GlassModule_dy*17/2." z="GlassModule_z0" begin="9" nmods="15"/>
- <line axis="x" x="GlassModule_dx/2." y="GlassModule_dy*19/2." z="GlassModule_z0" begin="8" nmods="16"/>
- <line axis="x" x="GlassModule_dx/2." y="GlassModule_dy*21/2." z="GlassModule_z0" begin="7" nmods="16"/>
- <line axis="x" x="GlassModule_dx/2." y="GlassModule_dy*23/2." z="GlassModule_z0" begin="6" nmods="17"/>
- <line axis="x" x="GlassModule_dx/2." y="GlassModule_dy*25/2." z="GlassModule_z0" begin="3" nmods="19"/>
- <line axis="x" x="GlassModule_dx/2." y="GlassModule_dy*27/2." z="GlassModule_z0" begin="0" nmods="22"/>
- <line axis="x" x="GlassModule_dx/2." y="GlassModule_dy*29/2." z="GlassModule_z0" begin="0" nmods="21"/>
- <line axis="x" x="GlassModule_dx/2." y="GlassModule_dy*31/2." z="GlassModule_z0" begin="0" nmods="20"/>
- <line axis="x" x="GlassModule_dx/2." y="GlassModule_dy*33/2." z="GlassModule_z0" begin="0" nmods="20"/>
- <line axis="x" x="GlassModule_dx/2." y="GlassModule_dy*35/2." z="GlassModule_z0" begin="0" nmods="19"/>
- <line axis="x" x="GlassModule_dx/2." y="GlassModule_dy*37/2." z="GlassModule_z0" begin="0" nmods="19"/>
- <line axis="x" x="GlassModule_dx/2." y="GlassModule_dy*39/2." z="GlassModule_z0" begin="0" nmods="17"/>
- <line axis="x" x="GlassModule_dx/2." y="GlassModule_dy*41/2." z="GlassModule_z0" begin="0" nmods="15"/>
- <line axis="x" x="GlassModule_dx/2." y="GlassModule_dy*43/2." z="GlassModule_z0" begin="0" nmods="14"/>
- <line axis="x" x="GlassModule_dx/2." y="GlassModule_dy*45/2." z="GlassModule_z0" begin="0" nmods="12"/>
- <line axis="x" x="GlassModule_dx/2." y="GlassModule_dy*47/2." z="GlassModule_z0" begin="0" nmods="10"/>
- <line axis="x" x="GlassModule_dx/2." y="GlassModule_dy*49/2." z="GlassModule_z0" begin="0" nmods="8"/>
- <line axis="x" x="GlassModule_dx/2." y="GlassModule_dy*51/2." z="GlassModule_z0" begin="0" nmods="7"/>
+ <line axis="x" x="GlassModule_dx/2." y="GlassModule_dy*1/2." z="GlassModule_z0" begin="13" nmods="10"/>
+ <line axis="x" x="GlassModule_dx/2." y="GlassModule_dy*3/2." z="GlassModule_z0" begin="13" nmods="10"/>
+ <line axis="x" x="GlassModule_dx/2." y="GlassModule_dy*5/2." z="GlassModule_z0" begin="13" nmods="10"/>
+ <line axis="x" x="GlassModule_dx/2." y="GlassModule_dy*7/2." z="GlassModule_z0" begin="12" nmods="10"/>
+ <line axis="x" x="GlassModule_dx/2." y="GlassModule_dy*9/2." z="GlassModule_z0" begin="12" nmods="10"/>
+ <line axis="x" x="GlassModule_dx/2." y="GlassModule_dy*11/2." z="GlassModule_z0" begin="12" nmods="10"/>
+ <line axis="x" x="GlassModule_dx/2." y="GlassModule_dy*13/2." z="GlassModule_z0" begin="11" nmods="11"/>
+ <line axis="x" x="GlassModule_dx/2." y="GlassModule_dy*15/2." z="GlassModule_z0" begin="10" nmods="11"/>
+ <line axis="x" x="GlassModule_dx/2." y="GlassModule_dy*17/2." z="GlassModule_z0" begin="9" nmods="12"/>
+ <line axis="x" x="GlassModule_dx/2." y="GlassModule_dy*19/2." z="GlassModule_z0" begin="8" nmods="13"/>
+ <line axis="x" x="GlassModule_dx/2." y="GlassModule_dy*21/2." z="GlassModule_z0" begin="7" nmods="13"/>
+ <line axis="x" x="GlassModule_dx/2." y="GlassModule_dy*23/2." z="GlassModule_z0" begin="6" nmods="14"/>
+ <line axis="x" x="GlassModule_dx/2." y="GlassModule_dy*25/2." z="GlassModule_z0" begin="3" nmods="16"/>
+ <line axis="x" x="GlassModule_dx/2." y="GlassModule_dy*27/2." z="GlassModule_z0" begin="0" nmods="18"/>
+ <line axis="x" x="GlassModule_dx/2." y="GlassModule_dy*29/2." z="GlassModule_z0" begin="0" nmods="18"/>
+ <line axis="x" x="GlassModule_dx/2." y="GlassModule_dy*31/2." z="GlassModule_z0" begin="0" nmods="16"/>
+ <line axis="x" x="GlassModule_dx/2." y="GlassModule_dy*33/2." z="GlassModule_z0" begin="0" nmods="16"/>
+ <line axis="x" x="GlassModule_dx/2." y="GlassModule_dy*35/2." z="GlassModule_z0" begin="0" nmods="14"/>
+ <line axis="x" x="GlassModule_dx/2." y="GlassModule_dy*37/2." z="GlassModule_z0" begin="0" nmods="13"/>
+ <line axis="x" x="GlassModule_dx/2." y="GlassModule_dy*39/2." z="GlassModule_z0" begin="0" nmods="11"/>
+ <line axis="x" x="GlassModule_dx/2." y="GlassModule_dy*41/2." z="GlassModule_z0" begin="0" nmods="10"/>
+ <line axis="x" x="GlassModule_dx/2." y="GlassModule_dy*43/2." z="GlassModule_z0" begin="0" nmods="7"/>
+ <line axis="x" x="GlassModule_dx/2." y="GlassModule_dy*45/2." z="GlassModule_z0" begin="0" nmods="3"/>
</lines>
</placements>
</detector>
diff --git a/compact/ci_ecal_scfi.xml b/compact/ci_ecal_scfi.xml
new file mode 100644
index 0000000000000000000000000000000000000000..a55596439d6146055952dd457bec2acc91eb1903
--- /dev/null
+++ b/compact/ci_ecal_scfi.xml
@@ -0,0 +1,59 @@
+<lccdd>
+ <define>
+ <constant name="EcalEndcapP_rmax" value="Solenoid_rmax"/>
+ <constant name="EcalEndcapP_FiberRadius" value="0.235*mm"/>
+ <constant name="EcalEndcapP_FiberOffset" value="0.5*mm"/>
+ <constant name="EcalEndcapP_FiberSpaceX" value="0.265*mm"/>
+ <constant name="EcalEndcapP_FiberSpaceY" value="0.425*mm"/>
+ </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="ScFiCalorimeter"
+ vis="InvisibleWithDaughters"
+ readout="EcalEndcapPHits">
+ <position x="0" y="0" z="EcalEndcapP_zmin + EcalEndcapP_length/2."/>
+ <dimensions rmin="EcalEndcapP_rmin" rmax="EcalEndcapP_rmax" length="EcalEndcapP_length"/>
+ <module sizex="25*mm" sizey="25*mm" sizez="170*mm" material="TungstenDens24" vis="EcalEndcapVis">
+ <fiber material="Polystyrene"
+ radius="EcalEndcapP_FiberRadius"
+ offset="EcalEndcapP_FiberOffset"
+ spacex="EcalEndcapP_FiberSpaceX"
+ spacey="EcalEndcapP_FiberSpaceY"/>
+ </module>
+ </detector>
+ </detectors>
+
+ <!-- Definition of the readout segmentation/definition -->
+ <readouts>
+ <readout name="EcalEndcapPHits">
+ <segmentation type="NoSegmentation"/>
+ <id>system:8,ring:8,module:20,fiber_x:8,fiber_y:8</id>
+ </readout>
+ </readouts>
+
+ <plugins>
+ </plugins>
+
+</lccdd>
+
diff --git a/compact/ecal.xml b/compact/ecal.xml
index 4f54d2058245118836b03b37f9164b98f1c9f326..0b743c8854bf7e2dda8286e3d3775597555bb3bb 100644
--- a/compact/ecal.xml
+++ b/compact/ecal.xml
@@ -8,13 +8,15 @@
<documentation level="10">
### Ecal configuration
</documentation>
- <include ref="ci_ecal.xml"/>
- <!--<include ref="compact/ci_ecal_shashlik.xml"/>-->
- <!--<include ref="compact/ce_ecal.xml"/>-->
-<!-- <include ref="ce_ecal_crystal_glass.xml"/>-->
+ <include ref="ci_ecal_scfi.xml"/>
+ <!--<include ref="ci_ecal.xml"/>-->
+
+ <!--<include ref="ce_ecal.xml"/>-->
+ <!--<include ref="ce_ecal_crystal_glass.xml"/>-->
<include ref="hybrid_ecal.xml"/>
- <!-- <include ref="compact/ecal_barrel.xml"/> -->
- <!-- <include ref="compact/ecal_barrel_hybrid.xml"/> -->
+
+ <!-- <include ref="ecal_barrel.xml"/> -->
+ <!-- <include ref="ecal_barrel_hybrid.xml"/> -->
<include ref="ecal_barrel_interlayers.xml"/>
</lccdd>
diff --git a/src/ScFiCalorimeter_geo.cpp b/src/ScFiCalorimeter_geo.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..5fc14cef776e270a52bb3b9e5c2fc5ede9756dc4
--- /dev/null
+++ b/src/ScFiCalorimeter_geo.cpp
@@ -0,0 +1,174 @@
+//==========================================================================
+// 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 shashlik 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;
+ for (int i = 0; i < nas; ++i) {
+ Assembly assembly(detName + Form("_ring%d", i + 1));
+ auto assemblyPV = env.placeVolume(assembly, Position{0., 0., 0.});
+ 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))));
+ }
+
+ 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);
+ }
+ }
+
+ return std::make_tuple(modVol, Position{sx, sy, sz});
+}
+
+DECLARE_DETELEMENT(ScFiCalorimeter, create_detector)
+