From 670ca618fa6b6558179985008effca996b26556a Mon Sep 17 00:00:00 2001
From: Dmitry Romanov <romanov@jlab.org>
Date: Fri, 13 Aug 2021 18:40:52 +0000
Subject: [PATCH] Resolve "Native parametrization of Negative Hybrid
Calorimeter"
---
CMakeLists.txt | 9 +-
compact/definitions.xml | 2 +-
compact/ecal.xml | 3 +-
compact/hybrid_ecal.xml | 73 +++++++++++
src/HybridCalorimeter_geo.cpp | 237 ++++++++++++++++++++++++++++++++++
5 files changed, 320 insertions(+), 4 deletions(-)
create mode 100644 compact/hybrid_ecal.xml
create mode 100644 src/HybridCalorimeter_geo.cpp
diff --git a/CMakeLists.txt b/CMakeLists.txt
index bfef7ec7..854bd0b2 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -6,10 +6,14 @@ PROJECT(athena
)
set(CMAKE_CXX_STANDARD 17)
-find_package( DD4hep REQUIRED COMPONENTS DDCore DDG4 )
+find_package( DD4hep REQUIRED COMPONENTS DDCore DDG4)
find_package(Acts REQUIRED COMPONENTS Core PluginIdentification PluginTGeo PluginDD4hep )
+find_package(fmt)
+#find_library(FMT_LIBRARY fmt)
+
+
#-----------------------------------------------------------------------------------
set(a_lib_name athena)
@@ -28,6 +32,7 @@ dd4hep_add_plugin(${a_lib_name} SOURCES
src/GaseousRICH_geo.cpp
src/GeometryHelpers.cpp
src/HomogeneousCalorimeter_geo.cpp
+ src/HybridCalorimeter_geo.cpp
src/MRich_geo.cpp
src/PolyhedraEndcapCalorimeter2_geo.cpp
src/ShashlikCalorimeter_geo.cpp
@@ -38,7 +43,7 @@ dd4hep_add_plugin(${a_lib_name} SOURCES
USES ActsCore ActsPluginDD4hep
)
target_link_libraries(${a_lib_name}
- PUBLIC DD4hep::DDCore DD4hep::DDRec
+ PUBLIC DD4hep::DDCore DD4hep::DDRec fmt::fmt
)
#-----------------------------------------------------------------------------------
diff --git a/compact/definitions.xml b/compact/definitions.xml
index 162f937e..9f96a84b 100644
--- a/compact/definitions.xml
+++ b/compact/definitions.xml
@@ -486,7 +486,7 @@ end of the solenoid coils.
EcalEndcapP_rmin and EcalEndcapN_rmin need to be set in sync with the forward and backward detectors
</documentation>
<constant name="EcalEndcapP_rmin" value="200.0*mm"/>
- <constant name="EcalEndcapN_rmin" value="300.0*mm"/>
+ <constant name="EcalEndcapN_rmin" value="100.0*mm"/>
<constant name="HcalEndcapP_rmin" value="EcalEndcapP_rmin"/>
<constant name="HcalEndcapN_rmin" value="EcalEndcapN_rmin"/>
diff --git a/compact/ecal.xml b/compact/ecal.xml
index de9ed8c6..4f54d205 100644
--- a/compact/ecal.xml
+++ b/compact/ecal.xml
@@ -11,7 +11,8 @@
<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="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_interlayers.xml"/>
diff --git a/compact/hybrid_ecal.xml b/compact/hybrid_ecal.xml
new file mode 100644
index 00000000..d53ced04
--- /dev/null
+++ b/compact/hybrid_ecal.xml
@@ -0,0 +1,73 @@
+<lccdd>
+ <define>
+ <comment>
+ Transition area.
+ The idea behind this parametrization is that:
+ one glass module with its wrap is always
+ a size of 4 crystal modules with its wraps.
+ Then the transition area (where glass meets crystals) has no gaps
+
+ +----------------+----------------+
+ | +----+ +----+ | +----------+ |
+ | | | | | | | | |
+ | +----+ +----+ | | | |
+ | +----+ +----+ | | | |
+ | | | | | | | | |
+ | +----+ +----+ | +----------+ |
+ +----------------+----------------+
+ crystal glass
+
+ This implies that:
+ GlassModule_wrap = 2*CrystalModule_wrap
+ GlassModule_sx = 2*CrystalModule_sx
+ GlassModule_sy = 2*CrystalModule_sy
+
+ </comment>
+
+ <constant name="CrystalModule_width" value="20.00*mm"/>
+ <constant name="CrystalModule_length" value="200.00*mm"/>
+ <constant name="CrystalModule_wrap" value="0.50*mm"/>
+ <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_wrap" value="2*CrystalModule_wrap"/>
+ <constant name="GlassModule_z0" value="0.0*cm"/>
+
+ <constant name="EcalEndcapN_rmax" value="80*cm"/>
+ <constant name="EcalEndcapN_rmin" value="5*cm"/>
+ <constant name="EcalEndcapN_thickness" value="GlassModule_length"/>
+ <constant name="EcalEndcapN_z0" value="-EcalEndcapN_zmin - EcalEndcapN_thickness/2"/>
+ <constant name="EcalEndcapNCrystal_rmax" value="40*cm"/>
+
+ <constant name="CrystalModule_distance" value="CrystalModule_width + CrystalModule_wrap"/>
+ <constant name="GlassModule_distance" value="GlassModule_width + GlassModule_wrap"/>
+ </define>
+
+ <display>
+ <vis name="HybridEcalOuterVis" alpha="0.5" r= "0.3" g="0.3" b="0.3" showDaughters="true" visible="true"/>
+ </display>
+
+ <detectors>
+
+ <documentation level="10">
+ #### Backwards Endcap EM Calorimeter
+
+ Backwards Endcap EM Calorimeter, placements generated by script
+ </documentation>
+ <detector id="ECalEndcapN_ID" name="EcalEndcapN" type="HybridCalorimeter" readout="EcalEndcapNHits">
+ <position x="0" y="0" z="EcalEndcapN_z0"/>
+ <rotation x="0" y="0" z="0"/>
+ </detector>
+ </detectors>
+ <readouts>
+ <comment>Effectively no segmentation, the segmentation is used to provide cell dimension info</comment>
+ <readout name="EcalEndcapNHits">
+ <segmentation type="MultiSegmentation" key="sector">
+ <segmentation name="CrystalSeg" key_value="1" type="CartesianGridXY" grid_size_x="CrystalModule_distance" grid_size_y="CrystalModule_distance"/>
+ <segmentation name="GlassSeg" key_value="2" type="CartesianGridXY" grid_size_x="GlassModule_distance" grid_size_y="GlassModule_distance"/>
+ </segmentation>
+ <id>system:8,sector:4,module:20,x:32:-16,y:-16</id>
+ </readout>
+ </readouts>
+</lccdd>
diff --git a/src/HybridCalorimeter_geo.cpp b/src/HybridCalorimeter_geo.cpp
new file mode 100644
index 00000000..58869213
--- /dev/null
+++ b/src/HybridCalorimeter_geo.cpp
@@ -0,0 +1,237 @@
+//==========================================================================
+// A general implementation for homogeneous calorimeter
+// it supports three types of placements
+// 1. Individual module placement with module dimensions and positions
+// 2. Array placement with module dimensions and numbers of row and columns
+// 3. Disk placement with module dimensions and (Rmin, Rmax), and (Phimin, Phimax)
+// 4. Lines placement with module dimensions and (mirrorx, mirrory)
+// (NOTE: anchor point is the 0th block of the line instead of line center)
+//--------------------------------------------------------------------------
+// Author: Chao Peng (ANL)
+// Date: 06/09/2021
+//==========================================================================
+
+#include "GeometryHelpers.h"
+#include "DD4hep/DetFactoryHelper.h"
+#include <XML/Helper.h>
+#include <iostream>
+#include <algorithm>
+#include <tuple>
+#include <math.h>
+#include <fmt/core.h>
+
+using namespace dd4hep;
+
+// main
+static Ref_t create_detector(Detector& desc, xml::Handle_t handle, SensitiveDetector sens)
+{
+
+ using namespace std;
+ using namespace fmt;
+
+ xml::DetElement detElem = handle;
+ std::string detName = detElem.nameStr();
+ int detID = detElem.id();
+ DetElement det(detName, detID);
+ sens.setType("calorimeter");
+
+ auto glass_material = desc.material("PbGlass");
+ auto crystal_material = desc.material("PbWO4");
+ auto air_material = desc.material("Air");
+
+ double ROut = desc.constantAsDouble("EcalEndcapN_rmax");
+ double RIn = desc.constantAsDouble("EcalEndcapN_rmin");
+ double SizeZ = desc.constantAsDouble("EcalEndcapN_thickness");
+ double glass_shift_z = desc.constantAsDouble("GlassModule_z0");
+ double Thickness = desc.constantAsDouble("EcalEndcapN_thickness");
+ double trans_radius = desc.constantAsDouble("EcalEndcapNCrystal_rmax");
+ double Glass_ShiftZ = desc.constantAsDouble("GlassModule_z0");
+ double Glass_Width = desc.constantAsDouble("GlassModule_width");
+ double Glass_Thickness = desc.constantAsDouble("GlassModule_length");
+ double Glass_Gap = desc.constantAsDouble("GlassModule_wrap");
+ double glass_distance = desc.constantAsDouble("GlassModule_distance");
+
+ double Crystal_Width = desc.constantAsDouble("CrystalModule_width");
+ double Crystal_Thickness = desc.constantAsDouble("CrystalModule_length");
+ double Crystal_Gap = desc.constantAsDouble("CrystalModule_wrap");
+ double crystal_distance = desc.constantAsDouble("CrystalModule_distance");
+ double Crystal_shift_z = desc.constantAsDouble("CrystalModule_z0");
+
+ // RIn and ROut will define outer tube embedding the calorimeter
+ // centers_rin/out define the maximum radius of module centers
+ // so that modules are not overlapping with mother tube volume
+ double hypotenuse = sqrt(0.5 * glass_distance * glass_distance);
+ double centers_rin = RIn + hypotenuse + 1*mm;
+ double centers_rout = ROut - hypotenuse - 1*mm;
+
+ // envelope
+
+ // Assembly assembly(detName);
+
+ Tube outer_tube(RIn, ROut, SizeZ / 2.0, 0., 360.0 * deg);
+ Volume ecal_vol("negative_ecal", outer_tube, air_material);
+ ecal_vol.setVisAttributes(desc.visAttributes("HybridEcalOuterVis"));
+
+ double Glass_OuterR = ROut - 1 * cm ;
+ double Glass_InnerR = trans_radius;
+ glass_shift_z = Thickness / 2. - Glass_Thickness / 2.;
+
+ // Geometry of modules
+ Box glass_box("glass_box", Glass_Width * 0.5, Glass_Width * 0.5, Glass_Thickness * 0.5);
+ Volume glass_module("glass_module", glass_box, glass_material);
+ glass_module.setVisAttributes(desc.visAttributes("EcalEndcapNModuleVis"));
+
+ Box crystal_box("crystal_box", Crystal_Width* 0.5, Crystal_Width * 0.5, Crystal_Thickness * 0.5);
+ Volume crystal_module("crystal_module", crystal_box, crystal_material);
+ crystal_module.setVisAttributes(desc.visAttributes("EcalEndcapNModuleVis"));
+
+ // GLASS
+ double diameter = 2 * Glass_OuterR;
+
+ // How many towers do we have per row/columnt?
+ // Add a gap + diameter as if we have N towers, we have N-1 gaps;
+ int towersInRow = std::ceil((diameter + Glass_Gap) / (Glass_Width + Glass_Gap));
+
+ // Is it odd or even number of towersInRow
+ double leftTowerPos, topTowerPos;
+ if(towersInRow%2) {
+ // |
+ // [ ][ ][ ][ ][ ]
+ // ^ |
+ int towersInHalfRow = std::ceil(towersInRow/2.0);
+ topTowerPos = leftTowerPos = -towersInHalfRow * (Glass_Width + Glass_Gap);
+
+ } else {
+ // |
+ // [ ][ ][ ][ ][ ][ ]
+ // ^ |
+ int towersInHalfRow = towersInRow/2;
+ topTowerPos = leftTowerPos = -(towersInHalfRow - 0.5) * (Glass_Width + Glass_Gap);
+ }
+
+ int moduleIndex = 0;
+
+// fmt::print("\nCE EMCAL GLASS SQUARE START\n");
+// fmt::print("Glass_Thickness = {} cm;\n", Glass_Thickness / cm);
+// fmt::print("Glass_Width = {} cm;\n", Glass_Width / cm);
+// fmt::print("Glass_Gap = {} cm;\n", Glass_Gap / cm);
+// fmt::print("Glass_InnerR = {} cm;\n", Glass_InnerR / cm);
+// fmt::print("Glass_OuterR = {} cm;\n", Glass_OuterR / cm);
+// fmt::print("Glass_PosZ = {} cm;\n", glass_shift_z / cm);
+// fmt::print("Towers in Row/Col = {} cm;\n", glass_shift_z / cm);
+// fmt::print("Top left tower pos = {:<10} {:<10} cm;\n", -leftTowerPos / cm, topTowerPos / cm);
+// fmt::print("#Towers info:\n");
+// fmt::print("#{:<5} {:<6} {:<3} {:<3} {:>10} {:>10} {}\n", "idx", "code", "col", "row", "x", "y", "name");
+
+
+ // We first do a "dry run", not really placing modules,
+ // but figuring out the ID scheme, number of modules, etc.
+ int glassModuleCount = 0;
+ int crystalModuleCount = 0;
+ int firstCrystRow = 1000000; // The first row, where crystals are started
+ int firstCrystCol = 1000000; // The fist column, where crystals are started
+ for(int rowIndex=0; rowIndex < towersInRow; rowIndex++) {
+ for(int colIndex=0; colIndex < towersInRow; colIndex++) {
+ double glass_x = leftTowerPos + colIndex * glass_distance;
+ double glass_y = topTowerPos + rowIndex * glass_distance;
+ double r = sqrt(glass_x * glass_x + glass_y * glass_y);
+
+ if (r < centers_rout && r > centers_rin) {
+ // we are placing something
+ if(r<trans_radius) {
+ // 4 crystal modules will be placed
+ crystalModuleCount+=4;
+
+ // Finding the first col and row where crystals start
+ // is the same algorithm as finding a minimum in array
+ if(colIndex<firstCrystCol) {
+ firstCrystCol = colIndex;
+ }
+ if(rowIndex<firstCrystRow) {
+ firstCrystRow = rowIndex;
+ }
+ }
+ else
+ {
+ // glass module will be places
+ glassModuleCount++;
+ }
+ }
+ }
+ }
+ // fmt::print("#Towers info:\n");
+ // fmt::print("#{:<5} {:<6} {:<3} {:<3} {:>10} {:>10} {}\n", "idx", "code", "col", "row", "x", "y", "name");
+ int glass_module_index = 0;
+ int cryst_module_index = 0;
+ for(int rowIndex=0; rowIndex < towersInRow; rowIndex++) {
+ for(int colIndex=0; colIndex < towersInRow; colIndex++) {
+ double glass_x = leftTowerPos + colIndex * (Glass_Width + Glass_Gap);
+ double glass_y = topTowerPos + rowIndex * (Glass_Width + Glass_Gap);
+ double r = sqrt(glass_x * glass_x + glass_y * glass_y);
+
+ if (r < centers_rout && r > centers_rin) {
+ int code = 1000 * rowIndex + colIndex;
+ std::string name = fmt::format("ce_EMCAL_glass_phys_{}", code);
+
+ if(r<trans_radius) {
+
+ // first crystal module
+ double crystal_x = glass_x - crystal_distance / 2;
+ double crystal_y = glass_y - crystal_distance / 2;
+ auto placement = ecal_vol.placeVolume(crystal_module, Position(crystal_x, crystal_y, Crystal_shift_z));
+ placement.addPhysVolID("sector", 1);
+ placement.addPhysVolID("module", cryst_module_index++);
+
+ // second crystal module
+ crystal_x = glass_x + crystal_distance / 2;
+ crystal_y = glass_y - crystal_distance / 2;
+ placement = ecal_vol.placeVolume(crystal_module, Position(crystal_x, crystal_y, Crystal_shift_z));
+ placement.addPhysVolID("sector", 1);
+ placement.addPhysVolID("module", cryst_module_index++);
+
+ // third crystal module
+ crystal_x = glass_x - crystal_distance / 2;
+ crystal_y = glass_y + crystal_distance / 2;
+ placement = ecal_vol.placeVolume(crystal_module, Position(crystal_x, crystal_y, Crystal_shift_z));
+ placement.addPhysVolID("sector", 1);
+ placement.addPhysVolID("module", cryst_module_index++);
+
+ // forth crystal module
+ crystal_x = glass_x + crystal_distance / 2;
+ crystal_y = glass_y + crystal_distance / 2;
+ placement = ecal_vol.placeVolume(crystal_module, Position(crystal_x, crystal_y, Crystal_shift_z));
+ placement.addPhysVolID("sector", 1);
+ placement.addPhysVolID("module", cryst_module_index++);
+ }
+ else
+ {
+ // glass module
+ auto placement = ecal_vol.placeVolume(glass_module, Position(glass_x, glass_y, glass_shift_z));
+ placement.addPhysVolID("sector", 2);
+ placement.addPhysVolID("module", glass_module_index++);
+ }
+
+
+ // fmt::print(" {:<5} {:<6} {:<3} {:<3} {:>10.4f} {:>10.4f} {}\n", towerIndex, code, colIndex, rowIndex, x / cm, y / cm, name);
+ glass_module_index++;
+ }
+ }
+ }
+// fmt::print("Total Glass modules: {}\n", towerIndex);
+// fmt::print("CE EMCAL GLASS END\n\n");
+
+ // detector position and rotation
+ auto pos = detElem.position();
+ auto rot = detElem.rotation();
+ Volume motherVol = desc.pickMotherVolume(det);
+ Transform3D tr(RotationZYX(rot.z(), rot.y(), rot.x()), Position(pos.x(), pos.y(), pos.z()));
+ PlacedVolume envPV = motherVol.placeVolume(ecal_vol, tr);
+ envPV.addPhysVolID("system", detID);
+ det.setPlacement(envPV);
+ return det;
+
+}
+
+//@}
+DECLARE_DETELEMENT(HybridCalorimeter, create_detector)
+
--
GitLab