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

implement Shashlik Calorimeter

parent cdebba0e
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 462 additions and 142 deletions
athena.xml
+ 0
− 6
View file @ 15115f46
...@@ -108,8 +108,6 @@ ...@@ -108,8 +108,6 @@
</display> </display>
<comment> Include the IP components first </comment> <comment> Include the IP components first </comment>
<include ref="ip6/forward_ion_beamline.xml"/>
<include ref="ip6/beampipe.xml"/>
<detectors> <detectors>
<detector id="VertexBarrelSubAssembly_ID" <detector id="VertexBarrelSubAssembly_ID"
...@@ -203,10 +201,6 @@ ...@@ -203,10 +201,6 @@
<include ref="compact/forward_trd.xml"/> <include ref="compact/forward_trd.xml"/>
<include ref="compact/gaseous_rich.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"/> <include ref="compact/mm_tracker_barrel.xml"/>
......
compact/ci_ecal.xml 0 → 100644
+ 76
− 0
View file @ 15115f46
<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>
compact/ci_ecal_shashlik.xml 0 → 100644
+ 76
− 0
View file @ 15115f46
<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>
compact/ecal.xml
+ 6
− 78
View file @ 15115f46
<lccdd> <lccdd>
<define> <comment>barrel ecal</comment>
<constant name="EcalEndcapP_rmax" value="Solenoid_rmax "/>
</define>
<limits>
</limits>
<regions>
</regions>
<!-- Common Generic visualization attributes -->
<comment>Common Generic visualization attributes</comment>
<display>
</display>
<!-- <include ref="ecal_barrel.xml"/> --> <!-- <include ref="ecal_barrel.xml"/> -->
<include ref="ecal_barrel_hybrid.xml"/> <include ref="ecal_barrel_hybrid.xml"/>
<comment>electron endcap ecal</comment>
<!--<include ref="ce_ecal.xml"/>--> <!--<include ref="ce_ecal.xml"/>-->
<include ref="ce_ecal_crystal_glass.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> <comment>hadron endcap ecal</comment>
</plugins> <include ref="ci_ecal.xml"/>
<!--<include ref="ci_ecal_shashlik.xml"/>-->
</lccdd> </lccdd>
src/GaseousRICH_geo.cpp
+ 1
− 1
View file @ 15115f46
...@@ -192,7 +192,7 @@ void build_sensors(Detector &desc, Volume &env, xml::Component plm, const Positi ...@@ -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.); double rotx = dd4hep::getAttrOrDefault<double>(sec, _Unicode(rotx), 0.);
// fill sensors to the piece // 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; int imod = 1;
for (auto &p : points) { for (auto &p : points) {
// transofrms are in a reversed order // transofrms are in a reversed order
......
src/GeometryHelpers.cpp
+ 97
− 37
View file @ 15115f46
#include "GeometryHelpers.h" #include "GeometryHelpers.h"
// some utility functions that can be shared // some utility functions that can be shared
namespace ref::utils { namespace athena::geo {
typedef ROOT::Math::XYPoint Point; 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 // 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) { if (pt.r() > rmax || pt.r() < rmin) {
return false; return false;
...@@ -27,33 +38,29 @@ inline bool in_ring(const Point &pt, double sx, double sy, double rmin, double r ...@@ -27,33 +38,29 @@ inline bool in_ring(const Point &pt, double sx, double sy, double rmin, double r
return true; 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 // 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 // std::cout << depth << "/" << max_depth << std::endl;
if (!in_ring(p, sx, sy, rmin, rmax, phmin, phmax) || overlap(p, sx, sy, res)) { // exceeds the maximum depth in searching or already placed
if ((depth > max_depth) || (already_placed(p, res, sx, sy))) {
return; 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 // continue search for a good placement or if no placement found yet
add_rectangle(Point(p.x() + sx, p.y()), res, sx, sy, rmin, rmax, phmin, phmax); if (in_ring || res.empty()) {
add_rectangle(Point(p.x() - sx, p.y()), res, sx, sy, rmin, rmax, phmin, phmax); // check adjacent squares
add_rectangle(Point(p.x(), p.y() + sy), res, sx, sy, rmin, rmax, phmin, phmax); 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); 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 // fill squares
...@@ -68,23 +75,76 @@ std::vector<Point> fillRectangles(Point ref, double sx, double sy, double rmin, ...@@ -68,23 +75,76 @@ std::vector<Point> fillRectangles(Point ref, double sx, double sy, double rmin,
// move to center // move to center
ref = ref - Point(int(ref.x()/sx)*sx, int(ref.y()/sy)*sy); 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) { std::vector<Point> res;
for (int ix = -n; ix < n; ++ix) { add_rectangle(ref, res, sx, sy, rmin, rmax, phmin, phmax, (int(rmax/sx) + 1)*(int(rmax/sy) + 1)*2);
for (int iy = -n; iy < n; ++iy) { return res;
Point pt(ref.x() + ix*sx, ref.y() + iy*sy); }
if (in_ring(pt, sx, sy, rmin, rmax, phmin, phmax)) {
return pt; // 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; std::vector<Point> res;
ref = find_seed(ref, int(rmax/sx) + 2, sx, sy, rmin, rmax, phmin, phmax); add_poly(ref, res, 6, lside, rmin, rmax, phmin, phmax, std::pow(int(rmax/lside) + 1, 2)*2);
add_rectangle(ref, res, sx, sy, rmin, rmax, phmin, phmax);
return res; return res;
} }
} // athena::geo
} // ref::utils
src/GeometryHelpers.h
+ 5
− 2
View file @ 15115f46
...@@ -3,7 +3,7 @@ ...@@ -3,7 +3,7 @@
#include "Math/Point2D.h" #include "Math/Point2D.h"
// some utility functions that can be shared // some utility functions that can be shared
namespace ref::utils { namespace athena::geo {
typedef ROOT::Math::XYPoint Point; typedef ROOT::Math::XYPoint Point;
...@@ -17,4 +17,7 @@ inline std::vector<Point> fillSquares(Point ref, double size, double rmin, doubl ...@@ -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); 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
src/HomogeneousCalorimeter_geo.cpp
+ 11
− 17
View file @ 15115f46
...@@ -16,10 +16,10 @@ ...@@ -16,10 +16,10 @@
#include <XML/Helper.h> #include <XML/Helper.h>
#include <iostream> #include <iostream>
#include <algorithm> #include <algorithm>
#include <tuple>
#include <math.h> #include <math.h>
using namespace dd4hep; using namespace dd4hep;
using namespace dd4hep::detail;
/** \addtogroup calorimeters Calorimeters /** \addtogroup calorimeters Calorimeters
*/ */
...@@ -170,13 +170,12 @@ static Ref_t create_detector(Detector& desc, xml::Handle_t handle, SensitiveDete ...@@ -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 // 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 mod = plm.child(_Unicode(module));
auto sx = mod.attr<double>(_Unicode(sizex)); auto sx = mod.attr<double>(_Unicode(sizex));
auto sy = mod.attr<double>(_Unicode(sizey)); auto sy = mod.attr<double>(_Unicode(sizey));
auto sz = mod.attr<double>(_Unicode(sizez)); auto sz = mod.attr<double>(_Unicode(sizez));
dim = Position{sx, sy, sz};
Box modShape(sx/2., sy/2., sz/2.); Box modShape(sx/2., sy/2., sz/2.);
auto modMat = desc.material(mod.attr<std::string>(_Unicode(material))); auto modMat = desc.material(mod.attr<std::string>(_Unicode(material)));
Volume modVol("module_vol", modShape, modMat); Volume modVol("module_vol", modShape, modMat);
...@@ -185,29 +184,27 @@ Volume build_module(Detector &desc, xml::Collection_t &plm, SensitiveDetector &s ...@@ -185,29 +184,27 @@ Volume build_module(Detector &desc, xml::Collection_t &plm, SensitiveDetector &s
// no wrapper // no wrapper
if (!plm.hasChild(_Unicode(wrapper))) { if (!plm.hasChild(_Unicode(wrapper))) {
return modVol; return std::make_tuple(modVol, Position{sx, sy, sz});
// build wrapper // build wrapper
} else { } else {
auto wrp = plm.child(_Unicode(wrapper)); auto wrp = plm.child(_Unicode(wrapper));
auto thickness = wrp.attr<double>(_Unicode(thickness)); auto thickness = wrp.attr<double>(_Unicode(thickness));
if (thickness == 0.) { if (thickness < 1e-12*mm) {
return modVol; return std::make_tuple(modVol, Position{sx, sy, sz});
} }
auto wrpMat = desc.material(wrp.attr<std::string>(_Unicode(material))); auto wrpMat = desc.material(wrp.attr<std::string>(_Unicode(material)));
Box wrpShape((sx + thickness)/2., (sy + thickness)/2., sz/2.); Box wrpShape((sx + thickness)/2., (sy + thickness)/2., sz/2.);
Volume wrpVol("wrapper_vol", wrpShape, wrpMat); Volume wrpVol("wrapper_vol", wrpShape, wrpMat);
wrpVol.placeVolume(modVol, Position(0., 0., 0.)); wrpVol.placeVolume(modVol, Position(0., 0., 0.));
wrpVol.setVisAttributes(desc.visAttributes(wrp.attr<std::string>(_Unicode(vis)))); wrpVol.setVisAttributes(desc.visAttributes(wrp.attr<std::string>(_Unicode(vis))));
dim = Position{sx + thickness, sy + thickness, sz}; return std::make_tuple(wrpVol, Position{sx + thickness, sy + thickness, sz});
return wrpVol;
} }
} }
// place modules, id must be provided // place modules, id must be provided
static void add_individuals(Detector& desc, Assembly &env, xml::Collection_t &plm, SensitiveDetector &sens, int sid) static void add_individuals(Detector& desc, Assembly &env, xml::Collection_t &plm, SensitiveDetector &sens, int sid)
{ {
Position modSize; auto [modVol, modSize] = build_module(desc, plm, sens);
auto modVol = build_module(desc, plm, sens, modSize);
int sector_id = dd4hep::getAttrOrDefault<int>(plm, _Unicode(sector), sid); int sector_id = dd4hep::getAttrOrDefault<int>(plm, _Unicode(sector), sid);
for (xml::Collection_t pl(plm, _Unicode(placement)); pl; ++pl) { 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 ...@@ -228,8 +225,7 @@ static void add_individuals(Detector& desc, Assembly &env, xml::Collection_t &pl
// place array of modules // place array of modules
static void add_array(Detector& desc, Assembly &env, xml::Collection_t &plm, SensitiveDetector &sens, int sid) static void add_array(Detector& desc, Assembly &env, xml::Collection_t &plm, SensitiveDetector &sens, int sid)
{ {
Position modSize; auto [modVol, modSize] = build_module(desc, plm, sens);
auto modVol = build_module(desc, plm, sens, modSize);
int sector_id = dd4hep::getAttrOrDefault<int>(plm, _Unicode(sector), sid); int sector_id = dd4hep::getAttrOrDefault<int>(plm, _Unicode(sector), sid);
int id_begin = dd4hep::getAttrOrDefault<int>(plm, _Unicode(id_begin), 1); int id_begin = dd4hep::getAttrOrDefault<int>(plm, _Unicode(id_begin), 1);
int nrow = plm.attr<int>(_Unicode(nrow)); int nrow = plm.attr<int>(_Unicode(nrow));
...@@ -266,8 +262,7 @@ static void add_array(Detector& desc, Assembly &env, xml::Collection_t &plm, Sen ...@@ -266,8 +262,7 @@ static void add_array(Detector& desc, Assembly &env, xml::Collection_t &plm, Sen
// place disk of modules // place disk of modules
static void add_disk(Detector& desc, Assembly &env, xml::Collection_t &plm, SensitiveDetector &sens, int sid) static void add_disk(Detector& desc, Assembly &env, xml::Collection_t &plm, SensitiveDetector &sens, int sid)
{ {
Position modSize; auto [modVol, modSize] = build_module(desc, plm, sens);
auto modVol = build_module(desc, plm, sens, modSize);
int sector_id = dd4hep::getAttrOrDefault<int>(plm, _Unicode(sector), sid); int sector_id = dd4hep::getAttrOrDefault<int>(plm, _Unicode(sector), sid);
int id_begin = dd4hep::getAttrOrDefault<int>(plm, _Unicode(id_begin), 1); int id_begin = dd4hep::getAttrOrDefault<int>(plm, _Unicode(id_begin), 1);
double rmin = plm.attr<double>(_Unicode(rmin)); double rmin = plm.attr<double>(_Unicode(rmin));
...@@ -275,7 +270,7 @@ static void add_disk(Detector& desc, Assembly &env, xml::Collection_t &plm, Sens ...@@ -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 phimin = dd4hep::getAttrOrDefault<double>(plm, _Unicode(phimin), 0.);
double phimax = dd4hep::getAttrOrDefault<double>(plm, _Unicode(phimax), 2.*M_PI); 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 // placement to mother
auto pos = get_xml_xyz(plm, _Unicode(position)); auto pos = get_xml_xyz(plm, _Unicode(position));
auto rot = get_xml_xyz(plm, _Unicode(rotation)); 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 ...@@ -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) // 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) static void add_lines(Detector& desc, Assembly &env, xml::Collection_t &plm, SensitiveDetector &sens, int sid)
{ {
Position modSize; auto [modVol, modSize] = build_module(desc, plm, sens);
auto modVol = build_module(desc, plm, sens, modSize);
int sector_id = dd4hep::getAttrOrDefault<int>(plm, _Unicode(sector), sid); int sector_id = dd4hep::getAttrOrDefault<int>(plm, _Unicode(sector), sid);
int id_begin = dd4hep::getAttrOrDefault<int>(plm, _Unicode(id_begin), 1); int id_begin = dd4hep::getAttrOrDefault<int>(plm, _Unicode(id_begin), 1);
bool mirrorx = dd4hep::getAttrOrDefault<bool>(plm, _Unicode(mirrorx), false); bool mirrorx = dd4hep::getAttrOrDefault<bool>(plm, _Unicode(mirrorx), false);
......
src/ShashlikCalorimeter.cpp 0 → 100644
+ 189
− 0
View file @ 15115f46
//==========================================================================
// 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)
src/ce_MRICH.cpp
+ 1
− 1
View file @ 15115f46
...@@ -82,7 +82,7 @@ void addModules(Volume &mother, xml::DetElement &detElem, Detector &desc, Sensit ...@@ -82,7 +82,7 @@ void addModules(Volume &mother, xml::DetElement &detElem, Detector &desc, Sensit
modVol.setSensitiveDetector(sens); modVol.setSensitiveDetector(sens);
// place modules in the sectors (disk) // 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 // determine module direction, always facing z = 0
double roty = dims.z() > 0. ? M_PI/2. : -M_PI/2.; double roty = dims.z() > 0. ? M_PI/2. : -M_PI/2.;
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment