Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • EIC/detectors/athena
  • zwzhao/athena
  • FernandoTA/athena
  • palspeic/athena
4 results
Show changes
Commits on Source (4)
Hide whitespace changes
Inline Side-by-side
Showing
with 278 additions and 218 deletions
.gitlab-ci.yml
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
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
......
athena.xml
View file @ 3f8b7ec3
......@@ -172,7 +172,7 @@
<documentation level="10">
### PID detectors
</documentation>
<!--include ref="compact/pid_config_acadia.xml" /-->
<include ref="compact/pid_config_acadia.xml" />
<documentation level="10">
## Central calorimetry
......
bin/generate_prim_file
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
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
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
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>
......
compact/hybrid_ecal.xml
View file @ 3f8b7ec3
......@@ -30,7 +30,7 @@
<constant name="CrystalModule_z0" value="10.*cm"/>
<constant name="GlassModule_width" value="2*CrystalModule_width"/>
<constant name="GlassModule_length" value="40.00*cm"/>
<constant name="GlassModule_length" value="55.00*cm"/>
<constant name="GlassModule_wrap" value="2*CrystalModule_wrap"/>
<constant name="GlassModule_z0" value="0.0*cm"/>
......
scripts/test_matscan.cxx deleted 100644 → 0
View file @ 6e18a867
R__LOAD_LIBRARY(libDDCore.so)
// R__LOAD_LIBRARY(libActsPluginDD4hep.so)
R__LOAD_LIBRARY(libDDG4.so)
R__LOAD_LIBRARY(libDDG4IO.so)
#include "DD4hep/Detector.h"
#include "DD4hep/DetElement.h"
#include "DD4hep/Objects.h"
#include "DD4hep/Detector.h"
#include "DDG4/Geant4Data.h"
#include "DDRec/CellIDPositionConverter.h"
#include "DDRec/SurfaceManager.h"
#include "DDRec/Surface.h"
#include "TCanvas.h"
#include "TChain.h"
#include "TGeoMedium.h"
#include "TGeoManager.h"
#include "DDRec/MaterialScan.h"
#include "DDRec/MaterialManager.h"
#include "DD4hep/Detector.h"
#include "DD4hep/Printout.h"
#include "fmt/core.h"
#include <iostream>
#include <fstream>
// #include "Acts/Geometry/TrackingGeometry.hpp"
// #include "Acts/Geometry/TrackingVolume.hpp"
// #include "Acts/Plugins/DD4hep/ConvertDD4hepDetector.hpp"
using namespace dd4hep;
using namespace dd4hep::rec;
void test_matscan(const char* compact = "athena.xml", TString face="z"){
dd4hep::Detector& detector = dd4hep::Detector::getInstance();
detector.fromCompact(compact);
MaterialScan matscan(detector);
fmt::print("\n");
fmt::print("All detector subsystem names:\n");
for(const auto& d : detector.detectors() ) {
fmt::print(" {}\n", d.first);
}
return;
TString det_list[14]={"TrackerBarrel_Inner","TrackerBarrel_Outer","TrackerEndcapN_Inner","TrackerEndcapN_Outer","TrackerEndcapP_Inner","TrackerEndcapP_Outer","TrackerSubAssembly_Inner","TrackerSubAssembly_Outer","VertexBarrel","VertexBarrelSubAssembly","VertexEndcapN","VertexEndcapP","VertexEndcapSubAssembly","cb_DIRC"};
for (int dd=0;dd<14;dd++){
TString detname = det_list[dd];
matscan.setDetector(detname.Data());
double x0=0,y0=0,z0=0,x1,y1,z1; // cm
double epsilon=1e-4; // (mm) default 1e-4: Materials with a thickness smaller than epsilon (default 1e-4=1mu
const char* fmt1 = "%7.2f %7.2f %7.2f %5d %-20s %3.0f %8.3f %8.4f %11.4f %11.4f %10.3f %8.2f %11.6f %11.6f %7.2f %7.2f %7.2f\n";
const char* fmt2 = "%7.2f %7.2f %7.2f %5d %-20s %3.0f %8.3f %8.4f %11.6g %11.6g %10.3f %8.2f %11.6f %11.6f %7.2f %7.2f %7.2f\n";
// x1 = 100; y1 = 100; z1 = 100;
// y1 = 100;
double a1,a2,a3;
for(a3=-100;a3<101;a3=a3+200){
TString fname = Form("/global/u2/s/shujie/eic/output/matscan/%s_%s%g.dat",detname.Data(),face.Data(),a3);
FILE * pFile;
pFile = fopen (fname,"w");
for(a1=-100;a1<100;a1=a1+1){
for(a2=-100;a2<100;a2=a2+1){
if (face=="x"){
y1=a1; z1=a2; x1=a3;
}
else if (face=="y"){
x1=a1; z1=a2; y1=a3;
}
else if (face=="z"){
x1=a1; y1=a2; z1=a3;
}
Vector3D p0(x0, y0, z0), p1(x1, y1, z1);
Vector3D end, direction;
direction = (p1-p0).unit();
const auto& placements = matscan.scan(x0,y0,z0,x1,y1,z1,epsilon);
// matscan.print(x0,y0,z0,x1,y1,z1,epsilon);
// return;
// Vector3D end, direction;
double sum_x0 = 0;
double sum_lambda = 0;
double path_length = 0, total_length = 0;
// TString fname = Form("/global/u2/s/shujie/eic_dir/matscan_%g_%g_%g.dat",x1,y1,z1);
for (unsigned i=0;i<placements.size();i++){
TGeoMaterial* mat=placements[i].first->GetMaterial();
double length = placements[i].second;
double nx0 = length / mat->GetRadLen();
double nLambda = length / mat->GetIntLen();
sum_x0 += nx0;
sum_lambda += nLambda;
path_length += length;
total_length += length;
end = p0 + total_length * direction;
const char* fmt = mat->GetRadLen() >= 1e5 ? fmt2 : fmt1;
// fprintf(pFile, "%d\n",i+1);
fprintf(pFile, fmt,x1, y1, z1, i+1, mat->GetName(), mat->GetZ(), mat->GetA(),
mat->GetDensity(), mat->GetRadLen(), mat->GetIntLen(),
length, path_length, sum_x0, sum_lambda, end[0], end[1], end[2]);
// ::printf(fmt, i+1, mat->GetName(), mat->GetZ(), mat->GetA(),
// mat->GetDensity(), mat->GetRadLen(), mat->GetIntLen(),
// length, path_length, sum_x0, sum_lambda, end[0], end[1], end[2]);
// // mat->Print();
}
}
cout<<detname<<" "<<x1<<","<<y1<<","<<z1<<endl;
// cout<<x1<<","<<y1<<","<<z1<<": "<<placements.size()<<" "<<sum_x0<<" "<<total_length<<endl;
}
fclose (pFile);
}
}
}
/*
<!-- <include ref="ecal_barrel_hybrid.xml"/> -->
<TGeoManager::CountLevels>: max level = 5, max placements = 2796
Error in <TGeoVoxelFinder::SortAll>: Volume B0Tracker: Cannot make slices on any axis
Error in <TGeoVoxelFinder::SortAll>: Volume ForwardRomanPot_Station_1: Cannot make slices on any axis
{ B0APF_BeamlineMagnet
B0PF_BeamlineMagnet
B0Tracker
B1APF_BeamlineMagnet
B1PF_BeamlineMagnet
B2PF_BeamlineMagnet
BPFR1_BeamlineMagnet
BackwardTOF
BarrelTOF
BeamPipe
EcalEndcapN
EcalEndcapP
ForwardOffMTracker
ForwardRomanPot_Station_1
ForwardTOF
ForwardTRD
GaseousRICH
HcalBarrel
HcalEndcapN
HcalEndcapP
Q1APF_BeamlineMagnet
Q1BPF_BeamlineMagnet
Q2PF_BeamlineMagnet
QPFC1_BeamlineMagnet
QPFC2_BeamlineMagnet
QPFC3_BeamlineMagnet
QPFC4_BeamlineMagnet
QPFR1_BeamlineMagnet
QPFR2_BeamlineMagnet
SolenoidCoilBarrel
SolenoidCoilEndcapN
SolenoidCoilEndcapP
TOFSubAssembly
TrackerBarrel_Inner
TrackerBarrel_Outer
TrackerEndcapN_Inner
TrackerEndcapN_Outer
TrackerEndcapP_Inner
TrackerEndcapP_Outer
TrackerSubAssembly_Inner
TrackerSubAssembly_Outer
VertexBarrel
VertexBarrelSubAssembly
VertexEndcapN
VertexEndcapP
VertexEndcapSubAssembly
cb_DIRC
ce_MRICH
ffi_ZDC_ECAL
ffi_ZDC_HCAL}
*/
src/ScFiCalorimeter_geo.cpp 0 → 100644
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)