Resolve "Native parametrization of Negative Hybrid Calorimeter"

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
   )
 
 #-----------------------------------------------------------------------------------

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"/>

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"/>

compact/hybrid_ecal.xml

+<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>

src/HybridCalorimeter_geo.cpp

+//==========================================================================
+//  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)
+