Implement module-theta readout (with per-layer merging)

Detector/DetFCCeeECalInclined/compact/FCCee_ECalBarrel_thetamodulemerged.xml

@@ -0,0 +1,143 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<lccdd xmlns:compact="http://www.lcsim.org/schemas/compact/1.0"
+       xmlns:xs="http://www.w3.org/2001/XMLSchema"
+       xs:noNamespaceSchemaLocation="http://www.lcsim.org/schemas/compact/1.0/compact.xsd">
+
+  <info name="FCCee_ECalBarrel"
+        title="Settings for FCCee Inclined ECal Barrel Calorimeter"
+        author="M.Aleksa,J.Faltova,A.Zaborowska, V. Volkl"
+        url="no"
+        status="development"
+        version="1.0">
+    <comment>
+      Settings for the inclined EM calorimeter.
+      The barrel is filled with liquid argon. Passive material includes lead in the middle and steal on the outside, glued together.
+      Passive plates are inclined by a certain angle from the radial direction.
+      In between of two passive plates there is a readout.
+      Space between the plate and readout is of trapezoidal shape and filled with liquid argon.
+      Definition of sizes, visualization settings, readout and longitudinal segmentation are specified. 
+    </comment>
+  </info>
+
+  <define>
+    <!-- Inclination angle of the lead plates -->
+    <constant name="InclinationAngle" value="50*degree"/>
+    <!-- thickness of active volume between two absorber plates at barrel Rmin, measured perpendicular to the readout plate -->
+    <constant name="LArGapThickness" value="1.239749*2*mm"/>
+
+    <!-- Air margin, thicknesses of cryostat and LAr bath -->
+    <constant name="AirMarginThickness" value="54*mm"/>    <!-- Space holder for air gap between cryostat vessels --> 
+
+    <constant name="CryoBarrelFrontWarm" value="10*mm"/>   <!-- Al solid corresponding to 0.11 X0 -->
+    <constant name="CryoBarrelFrontCold" value="3.8*mm"/>  <!-- Al solid equivalent of 0.043 X0 sandwich CFRP -->
+    <constant name="CryoBarrelFront" value="CryoBarrelFrontWarm+CryoBarrelFrontCold"/>
+
+    <constant name="CryoBarrelBackCold" value="30*mm"/>   <!-- Al solid corresponding to 0.34 X0 -->
+    <constant name="CryoBarrelBackWarm" value="2.7*mm"/>  <!-- Al solid equivalent of 0.03 X0 sandwich CFRP -->
+    <constant name="SolenoidBarrel" value="70*mm"/>    <!-- Al solenoid with thickness of 0.8 X0 -->
+    <constant name="CryoBarrelBack" value="CryoBarrelBackWarm+SolenoidBarrel+CryoBarrelBackCold"/>
+
+    <constant name="CryoBarrelSideWarm" value="30*mm"/>
+    <constant name="CryoBarrelSideCold" value="3.8*mm"/>
+    <constant name="CryoBarrelSide" value="CryoBarrelSideWarm+CryoBarrelSideCold"/>
+
+    <constant name="LArBathThicknessFront" value="5*mm"/>
+    <constant name="LArBathThicknessBack" value="40*mm"/>
+
+    <!-- air margin around calorimeter -->
+    <constant name="BarCryoECal_rmin" value="BarECal_rmin+AirMarginThickness"/>
+    <constant name="BarCryoECal_rmax" value="BarECal_rmax-AirMarginThickness"/>
+    <constant name="BarCryoECal_dz" value="BarECal_dz"/>
+    <!-- calorimeter active volume -->
+    <constant name="EMBarrel_rmin" value="BarCryoECal_rmin+CryoBarrelFront+LArBathThicknessFront"/>
+    <constant name="EMBarrel_rmax" value="BarCryoECal_rmax-CryoBarrelBack-LArBathThicknessBack"/>
+    <constant name="EMBarrel_dz" value="BarECal_dz-CryoBarrelSide"/>
+    <!-- thickness of active volume between two absorber plates at EMBarrel_rmin, measuring perpendicular to the readout plate -->
+    <constant name="LAr_thickness" value="LArGapThickness"/>
+    <!-- passive layer consists of lead in the middle and steel on the outside, glued -->
+    <!-- When employing trapezoidal planes Pb_thickness corresponds to the minimum thickness, i.e at the front of the calo -->
+    <constant name="Pb_thickness" value="1.80*mm"/>
+    <constant name="planeLength" value="-EMBarrel_rmin*cos(InclinationAngle) + sqrt(EMBarrel_rmax*EMBarrel_rmax - EMBarrel_rmin*EMBarrel_rmin*sin(InclinationAngle)*sin(InclinationAngle))"/>
+    <constant name="ECalBarrelNumPlanes" value="1545"/>
+    <constant name="ECalBarrelNumLayers" value="12"/>
+    <constant name="phi" value="asin(planeLength / EMBarrel_rmax * sin(InclinationAngle))"/>
+    <!-- use a different value for Pb_thickness_max when employing trapezoidal planes -->
+    <!-- approximate constant sampling fraction: make the absorber grow linearly with the radius,
+         taking into account the angular projection effect -->
+    <!-- <constant name="Pb_thickness_max" value="1.3 * Pb_thickness * EMBarrel_rmax/EMBarrel_rmin *
+         cos(InclinationAngle - phi) / cos(InclinationAngle)" />-->
+    <constant name="Pb_thickness_max" value="Pb_thickness" />
+    <!-- total amount of steel in one passive plate: it is divided for the outside layer on top and bottom -->
+    <constant name="Steel_thickness" value="0.1*mm"/>
+    <!-- total amount of glue in one passive plate: it is divided for the outside layer on top and bottom -->
+    <constant name="Glue_thickness" value="0.1*mm"/>
+    <!-- readout in between two absorber plates -->
+    <constant name="readout_thickness" value="1.2*mm"/>
+  </define>
+
+  <display>
+    <vis name="ecal_envelope" r="0.1" g="0.2" b="0.6" alpha="1" showDaughers="false" visible="true" />
+  </display>
+
+  <readouts>
+    <!-- readout for the simulation, with the baseline merging: 2x along the module direction in each layer; 4x along theta in each layer except layer 1 -->
+    <!-- the lists mergedCells_Theta and mergedModules define the number of cells to group together in the theta and module direction as a function of the layer -->
+    <readout name="ECalBarrelModuleThetaMerged">
+      <segmentation type="FCCSWGridModuleThetaMerged" nModules="1545" mergedCells_Theta="4 1 4 4 4 4 4 4 4 4 4 4" mergedModules="2 2 2 2 2 2 2 2 2 2 2 2" grid_size_theta="0.009817477/4" offset_theta="0.5902785"/>
+        <id>system:4,cryo:1,type:3,subtype:3,layer:8,module:11,theta:10</id>
+    </readout>
+
+    <!-- example of adding a second readout for the reconstruction, to compare the two -->
+    <readout name="ECalBarrelModuleThetaMerged2">
+      <segmentation type="FCCSWGridModuleThetaMerged" nModules="1545" mergedCells_Theta="2 4 2 1 2 1 2 2 1 1 1 2" mergedModules="2 1 1 2 2 1 1 1 2 2 1 1" grid_size_theta="0.009817477/4" offset_theta="0.5902785"/>
+        <id>system:4,cryo:1,type:3,subtype:3,layer:8,module:11,theta:10</id>
+    </readout>
+  </readouts>
+
+  <detectors>
+    <detector id="BarECal_id" name="ECalBarrel" type="EmCaloBarrelInclined" readout="ECalBarrelModuleThetaMerged">
+      <type_flags type=" DetType_CALORIMETER + DetType_ELECTROMAGNETIC + DetType_BARREL"/>
+      <sensitive type="SimpleCalorimeterSD"/>
+      <dimensions rmin="BarCryoECal_rmin" rmax="BarCryoECal_rmax" dz="BarCryoECal_dz" vis="ecal_envelope"/>
+      <cryostat name="ECAL_Cryo">
+        <material name="Aluminum"/>
+	<dimensions rmin1="BarCryoECal_rmin" rmin2="BarCryoECal_rmin+CryoBarrelFront" rmax1="BarCryoECal_rmax-CryoBarrelBack" rmax2="BarCryoECal_rmax" dz="BarCryoECal_dz"/>
+        <front sensitive="false"/> <!-- inner wall of the cryostat -->
+        <side sensitive="false"/> <!-- both sides of the cryostat -->
+        <back sensitive="false"/> <!-- outer wall of the cryostat -->
+      </cryostat>
+      <calorimeter name="EM_barrel">
+        <!-- offset defines the numbering of the modules: module==0 for phi=0 direction -->
+        <dimensions rmin="EMBarrel_rmin" rmax="EMBarrel_rmax" dz="EMBarrel_dz" offset="-InclinationAngle"/>
+        <active thickness="LAr_thickness">
+          <material name="LAr"/>
+          <!-- overlap offset is a specific feature of the construction; do not change! -->
+          <!-- one volume for a gap on both side of the readout) --> 
+          <overlap offset="0.5"/>
+        </active>
+        <passive>
+          <rotation angle="InclinationAngle"/>  <!-- inclination angle -->
+          <inner thickness="Pb_thickness" sensitive="false">
+            <material name="Lead"/>
+          </inner>
+          <innerMax thickness="Pb_thickness_max" sensitive="false">
+            <material name="Lead"/>
+          </innerMax>
+          <glue thickness="Glue_thickness" sensitive="false">
+            <material name="lArCaloGlue"/>
+          </glue>
+          <outer thickness="Steel_thickness" sensitive="false">
+            <material name="lArCaloSteel"/>
+          </outer>
+        </passive>
+        <readout thickness="readout_thickness" sensitive="false">
+          <material name="PCB"/>
+        </readout>
+        <layers> <!-- pcb electrode segmentation in the radial direction -->
+           <layer thickness="1.5*cm" repeat="1"/>
+           <layer thickness="3.5*cm" repeat="11"/>
+        </layers>
+      </calorimeter>
+    </detector>
+  </detectors>
+</lccdd>

Detector/DetFCCeeIDEA-LAr/compact/FCCee_DectMaster_thetamodulemerged.xml

@@ -0,0 +1,45 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<lccdd xmlns:compact="http://www.lcsim.org/schemas/compact/1.0"
+  xmlns:xs="http://www.w3.org/2001/XMLSchema"
+  xs:noNamespaceSchemaLocation="http://www.lcsim.org/schemas/compact/1.0/compact.xsd">
+
+  <info name="FCCee-IDEA-LAr Master"
+    title="FCCee-IDEA-LAr Master"
+    author="Valentin Volkl"
+    url="no"
+    status="development"
+    version="1.0">
+    <comment>
+      Master compact file describing the latest developments of the FCCee IDEA detector concept with a LAr calorimeter.
+    </comment>
+  </info>
+
+  <include ref="${DD4hepINSTALL}/DDDetectors/compact/detector_types.xml" />
+
+  <includes>
+    <gdmlFile  ref="../../DetCommon/compact/elements.xml"/>
+    <gdmlFile  ref="../../DetCommon/compact/materials.xml"/>
+  </includes>
+
+  <define>
+    <constant name="world_size" value="25*m"/>
+    <constant name="world_x" value="world_size"/>
+    <constant name="world_y" value="world_size"/>
+    <constant name="world_z" value="world_size"/>
+  </define>
+
+  <include ref="./FCCee_DectDimensions.xml" />
+
+  <include ref="../../DetFCCeeIDEA/compact/Beampipe.xml"/>
+  <include ref="../../DetFCCeeIDEA/compact/BeamInstrumentation.xml"/>
+  <include ref="../../DetFCCeeIDEA/compact/LumiCal.xml"/>
+  <include ref="../../DetFCCeeIDEA/compact/HOMAbsorber.xml"/> 
+  <include ref="../../DetFCCeeCLD/compact/FCCee_o2_v02/Vertex.xml"/>
+  <include ref="../../DetFCCeeIDEA/compact/SimplifiedDriftChamber.xml"/>
+  <include ref="../../DetFCCeeECalInclined/compact/FCCee_ECalBarrel_thetamodulemerged.xml"/>
+  <include ref="../../DetFCCeeHCalTile/compact/FCCee_HCalBarrel_TileCal.xml"/>
+  <include ref="../../DetFCCeeCalDiscs/compact/FCCee_EcalEndcaps_coneCryo.xml"/>
+  <include ref="../../DetFCCeeHCalTile/compact/FCCee_HCalEndcaps_ThreeParts_TileCal.xml"/>
+  <include ref="MuonTagger.xml"/>
+
+</lccdd>

Detector/DetFCChhECalInclined/src/ECalBarrelInclined_geo.cpp

@@ -66,6 +66,25 @@ static dd4hep::detail::Ref_t createECalBarrelInclined(dd4hep::Detector& aLcdd,
     layersTotalHeight += layer.repeat() * layer.thickness();
   }
   lLog << MSG::DEBUG << "Number of layers: " << numLayers << " total thickness " << layersTotalHeight << endmsg;
+  // The following code checks if the xml geometry file contains a constant defining
+  // the number of layers the barrel. In that case, it makes the program abort
+  // if the number of planes in the xml is different from the one calculated from
+  // the geometry. This is because the number of layers is needed
+  // in other parts of the code (the readout for the FCC-ee ECAL with
+  // inclined modules).
+  int nLayers = -1;
+  try {
+    nLayers = aLcdd.constant<int>("ECalBarrelNumLayers");
+  }
+  catch(...) {
+    ;
+  }
+  if (nLayers > 0 && nLayers != numLayers) {
+    lLog << MSG::ERROR << "Incorrect number of layers (ECalBarrelNumLayers) in xml file!" << endmsg;
+    // todo: incidentSvc->fireIncident(Incident("ECalConstruction", "GeometryFailure"));
+    // make the code crash (incidentSvc does not work)
+     assert(nLayers == numLayers);
+  }
 
   dd4hep::xml::DetElement readout = calo.child(_Unicode(readout));
   std::string readoutMaterial = readout.materialStr();
@@ -185,6 +204,29 @@ static dd4hep::detail::Ref_t createECalBarrelInclined(dd4hep::Detector& aLcdd,
        << " rotation angle = " << angle << endmsg;
   uint numPlanes =
       round(M_PI / asin((passiveThickness + activeThickness + readoutThickness) / (2. * caloDim.rmin() * cos(angle))));
+
+  // The following code checks if the xml geometry file contains a constant defining
+  // the number of planes in the barrel. In that case, it makes the program abort
+  // if the number of planes in the xml is different from the one calculated from
+  // the geometry. This is because the number of plane information (retrieved from the
+  // xml) is used in other parts of the code (the readout for the FCC-ee ECAL with
+  // inclined modules). In principle the code above should be refactored so that the number
+  // of planes is one of the inputs of the calculation and other geometrical parameters
+  // are adjusted accordingly. This is left for the future, and we use the workaround
+  // below to enforce for the time being that the number of planes is "correct"
+  int nModules = -1;
+  try {
+    nModules = aLcdd.constant<int>("ECalBarrelNumPlanes");
+  }
+  catch(...) {
+    ;
+  }
+  if (nModules > 0 && nModules != numPlanes) {
+    lLog << MSG::ERROR << "Incorrect number of planes (ECalBarrelNumPlanes) in xml file!" << endmsg;
+    // todo: incidentSvc->fireIncident(Incident("ECalConstruction", "GeometryFailure"));
+    // make the code crash (incidentSvc does not work)
+     assert(nModules == numPlanes);
+  }
   double dPhi = 2. * M_PI / numPlanes;
   lLog << MSG::INFO << "number of passive plates = " << numPlanes << " azim. angle difference =  " << dPhi << endmsg;
   lLog << MSG::INFO << " distance at inner radius (cm) = " << 2. * M_PI * caloDim.rmin() / numPlanes << "\n"

Detector/DetSegmentation/include/DetSegmentation/FCCSWGridModuleThetaMerged.h

@@ -0,0 +1,117 @@
+#ifndef DETSEGMENTATION_FCCSWGRIDMODULETHETAMERGED_H
+#define DETSEGMENTATION_FCCSWGRIDMODULETHETAMERGED_H
+
+// FCCSW
+#include "DetSegmentation/GridTheta.h"
+#include "DD4hep/VolumeManager.h"
+
+/** FCCSWGridModuleThetaMerged Detector/DetSegmentation/DetSegmentation/FCCSWGridModuleThetaMerged.h FCCSWGridModuleThetaMerged.h
+ *
+ *  Segmentation in theta and module.
+ *  Based on GridTheta, merges modules and theta cells based on layer number
+ *
+ */
+
+namespace dd4hep {
+namespace DDSegmentation {
+class FCCSWGridModuleThetaMerged : public GridTheta {
+public:
+  /// default constructor using an arbitrary type
+  FCCSWGridModuleThetaMerged(const std::string& aCellEncoding);
+  /// Default constructor used by derived classes passing an existing decoder
+  FCCSWGridModuleThetaMerged(const BitFieldCoder* decoder);
+
+  /// destructor
+  virtual ~FCCSWGridModuleThetaMerged() = default;
+
+  /// read n(modules) from detector metadata
+  void GetNModulesFromGeom();
+
+  /// read n(layers) from detector metadata
+  void GetNLayersFromGeom();
+
+  /**  Determine the local position based on the cell ID.
+   *   @param[in] aCellId ID of a cell.
+   *   return Position (relative to R, phi of Geant4 volume it belongs to, scaled for R=1).
+   */
+  virtual Vector3D position(const CellID& aCellID) const;
+  /**  Determine the cell ID based on the position.
+   *   @param[in] aLocalPosition (not used).
+   *   @param[in] aGlobalPosition
+   *   @param[in] aVolumeId ID of the Geant4 volume
+   *   return Cell ID.
+   */
+  virtual CellID cellID(const Vector3D& aLocalPosition, const Vector3D& aGlobalPosition,
+                        const VolumeID& aVolumeID) const;
+  /**  Determine the azimuthal angle (relative to the G4 volume) based on the cell ID.
+   *   @param[in] aCellId ID of a cell.
+   *   return Phi.
+   */
+  double phi(const CellID& aCellID) const;
+  /**  Determine the polar angle (relative to the G4 volume) based on the cell ID.
+   *   @param[in] aCellId ID of a cell.
+   *   return Theta.
+   */
+  double theta(const CellID& aCellID) const;
+  /**  Determine the radius based on the cell ID.
+   *   @param[in] aCellId ID of a cell.
+   *   return Radius.
+   */
+  // double radius(const CellID& aCellID) const;
+  /**  Get the number of merged cells in theta for given layer
+   *   @param[in] layer
+   *   return The number of merged cells in theta
+   */
+  inline int mergedThetaCells(const int layer) const {
+    if (layer<m_mergedCellsTheta.size())
+      return m_mergedCellsTheta[layer];
+    else
+      return 1;
+  }
+  /**  Get the number of merged modules (inclined in phi)
+   *   @param[in] layer
+   *   return The number of merged modules
+   */
+  inline int mergedModules(const int layer) const {
+    if (layer<m_mergedModules.size())
+      return m_mergedModules[layer];
+    else
+      return 1;
+  }
+  /**  Get the total number of modules of detector
+   *   return The number of modules (as it was set by the user in the xml file..)
+   */
+  inline int nModules() const { return m_nModules; }
+  /**  Get the number of layers
+   *   return The number of layers
+   */
+  inline int nLayers() const { return m_nLayers; }
+  /**  Get the field name used for the layer
+   *   return The field name for the layer.
+   */
+  inline const std::string& fieldNameLayer() const { return m_layerID; }
+  /**  Get the field name used for the module number
+   *   return The field name for the module.
+   */
+  inline const std::string& fieldNameModule() const { return m_moduleID; }
+
+protected:
+  /// the field name used for layer
+  std::string m_layerID;
+  /// the field name used for the read-out module (can differ from module due to merging)
+  std::string m_moduleID;
+  /// vector of number of cells to be merged along theta for each layer
+  std::vector<int> m_mergedCellsTheta;
+  /// vector of number of modules to be merged for each layer
+  std::vector<int> m_mergedModules;
+
+  /// number of modules (or, equivalently, the deltaPhi between adjacent modules)
+  int m_nModules;
+
+  /// number of layers (from the geometry)
+  int m_nLayers;
+
+};
+}
+}
+#endif /* DETSEGMENTATION_FCCSWGRIDMODULETHETAMERGED_H */