Skip to content
Snippets Groups Projects
Commit 3f8b7ec3 authored by Chao Peng's avatar Chao Peng
Browse files

Resolve "Implement WScFi for the hadron endcap"

parent 15b7821c
Branches
Tags
1 merge request!118Resolve "Implement WScFi for the hadron endcap"
Hide whitespace changes
Inline Side-by-side
.gitlab-ci.yml
+ 2
− 0
View file @ 3f8b7ec3
......@@ -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
......
CMakeLists.txt
+ 1
− 0
View file @ 3f8b7ec3
......@@ -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
......
bin/generate_prim_file
+ 9
− 1
View file @ 3f8b7ec3
......@@ -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='')
......
compact/ce_ecal_crystal_glass.xml
+ 23
− 26
View file @ 3f8b7ec3
......@@ -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>
......
compact/ci_ecal_scfi.xml 0 → 100644
+ 59
− 0
View file @ 3f8b7ec3
<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>
compact/ecal.xml
+ 8
− 6
View file @ 3f8b7ec3
......@@ -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>
......
src/ScFiCalorimeter_geo.cpp 0 → 100644
+ 174
− 0
View file @ 3f8b7ec3
//==========================================================================
// 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)
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment