From 9542241ac18dcf56636361318d2dae088f5e7c70 Mon Sep 17 00:00:00 2001
From: Giovanni Marchiori <giovanni.marchiori@cern.ch>
Date: Thu, 21 Sep 2023 23:22:10 +0200
Subject: [PATCH] add number of layers info (currently only for printout of
 merging config but could be used to enforce that length of vector passed via
 xml to configure merging per layer is correct

---
 .../FCCee_ECalBarrel_thetamodulemerged.xml    |  1 +
 .../src/ECalBarrelInclined_geo.cpp            | 21 +++++++++++-
 .../FCCSWGridModuleThetaMerged.h              | 34 ++++++++++++++-----
 .../FCCSWGridModuleThetaMergedHandle.h        | 15 ++++----
 .../src/FCCSWGridModuleThetaMerged.cpp        | 18 ++++++++--
 5 files changed, 68 insertions(+), 21 deletions(-)

diff --git a/Detector/DetFCCeeECalInclined/compact/FCCee_ECalBarrel_thetamodulemerged.xml b/Detector/DetFCCeeECalInclined/compact/FCCee_ECalBarrel_thetamodulemerged.xml
index e9401b96..c9a8f312 100644
--- a/Detector/DetFCCeeECalInclined/compact/FCCee_ECalBarrel_thetamodulemerged.xml
+++ b/Detector/DetFCCeeECalInclined/compact/FCCee_ECalBarrel_thetamodulemerged.xml
@@ -59,6 +59,7 @@
     <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,
diff --git a/Detector/DetFCChhECalInclined/src/ECalBarrelInclined_geo.cpp b/Detector/DetFCChhECalInclined/src/ECalBarrelInclined_geo.cpp
index 72a7bb0c..31abfcfa 100644
--- a/Detector/DetFCChhECalInclined/src/ECalBarrelInclined_geo.cpp
+++ b/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();
@@ -203,7 +222,7 @@ static dd4hep::detail::Ref_t createECalBarrelInclined(dd4hep::Detector& aLcdd,
     ;
   }
   if (nModules > 0 && nModules != numPlanes) {
-    lLog << MSG::ERROR << "Incorrect number of planes in xml file!" << endmsg;
+    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);
diff --git a/Detector/DetSegmentation/include/DetSegmentation/FCCSWGridModuleThetaMerged.h b/Detector/DetSegmentation/include/DetSegmentation/FCCSWGridModuleThetaMerged.h
index 18e0a83f..c549eabe 100644
--- a/Detector/DetSegmentation/include/DetSegmentation/FCCSWGridModuleThetaMerged.h
+++ b/Detector/DetSegmentation/include/DetSegmentation/FCCSWGridModuleThetaMerged.h
@@ -24,8 +24,11 @@ class FCCSWGridModuleThetaMerged : public GridTheta {
   /// destructor
   virtual ~FCCSWGridModuleThetaMerged() = default;
 
-  /// read nmodules from detector metadata
+  /// 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.
@@ -59,20 +62,30 @@ class FCCSWGridModuleThetaMerged : public GridTheta {
    *   @param[in] layer
    *   return The number of merged cells in theta
    */
-  inline int mergedThetaCells(const int layer) const { return m_mergedCellsTheta[layer]; }
+  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 cells in theta
+   *   return The number of merged modules
    */
-  inline int mergedModules(const int layer) const { return m_mergedModules[layer]; }
-  /**  Get the number of 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; }
-  /**  Set the number of modules
-   *   return The number of modules (as it was set by the user in the xml file..)
+  /**  Get the number of layers
+   *   return The number of layers
    */
-  inline void setNModules(const int nModules) { m_nModules = nModules; }
+  inline int nLayers() const { return m_nLayers; }
   /**  Get the field name used for the layer
    *   return The field name for the layer.
    */
@@ -92,9 +105,12 @@ class FCCSWGridModuleThetaMerged : public GridTheta {
   /// vector of number of modules to be merged for each layer
   std::vector<int> m_mergedModules;
 
-  /// to be seen how to retrieve this initialization step from the geometry
   /// number of modules (or, equivalently, the deltaPhi between adjacent modules)
   int m_nModules;
+
+  /// number of layers (from the geometry)
+  int m_nLayers;
+  
 };
 }
 }
diff --git a/Detector/DetSegmentation/include/DetSegmentation/FCCSWGridModuleThetaMergedHandle.h b/Detector/DetSegmentation/include/DetSegmentation/FCCSWGridModuleThetaMergedHandle.h
index 0d4f561f..bf6b6bcb 100644
--- a/Detector/DetSegmentation/include/DetSegmentation/FCCSWGridModuleThetaMergedHandle.h
+++ b/Detector/DetSegmentation/include/DetSegmentation/FCCSWGridModuleThetaMergedHandle.h
@@ -70,7 +70,6 @@ class FCCSWGridModuleThetaMerged : public FCCSWGridModuleThetaMergedHandle {
     return access()->implementation->cellID(local, global, volID);
   }
 
-
   /// access the grid size in theta
   inline double gridSizeTheta() const { return access()->implementation->gridSizeTheta(); }
 
@@ -83,20 +82,18 @@ class FCCSWGridModuleThetaMerged : public FCCSWGridModuleThetaMergedHandle {
   /// access the number of modules
   inline int nModules() const { return access()->implementation->nModules(); }
 
+  /// access the number of layers
+  inline int nLayers() const { return access()->implementation->nLayers(); }
+
   /// access the number of merged modules for given layer
   inline int mergedModules(const int layer) const { return access()->implementation->mergedModules(layer); }
 
-
   /// set the coordinate offset in theta
   inline void setOffsetTheta(double offset) const { access()->implementation->setOffsetTheta(offset); }
 
   /// set the grid size in theta
   inline void setGridSizeTheta(double cellSize) const { access()->implementation->setGridSizeTheta(cellSize); }
 
-  /// set the number of modules
-  inline void setNModules(int nModules) const { access()->implementation->setNModules(nModules); }
-
-
   /// access the field name used for theta
   inline const std::string& fieldNameTheta() const { return access()->implementation->fieldNameTheta(); }
 
@@ -108,12 +105,12 @@ class FCCSWGridModuleThetaMerged : public FCCSWGridModuleThetaMergedHandle {
 
 
   /** \brief Returns a std::vector<double> of the cellDimensions of the given cell ID
-      in natural order of dimensions (dPhi, dTheta)
+      in natural order of dimensions (nModules, dTheta)
 
       Returns a std::vector of the cellDimensions of the given cell ID
-      \param cellID is ignored as all cells have the same dimension
+      \param cellID
       \return std::vector<double> size 2:
-      -# size in phi
+      -# size in module
       -# size in theta
   */
   inline std::vector<double> cellDimensions(const CellID& /*id*/) const {
diff --git a/Detector/DetSegmentation/src/FCCSWGridModuleThetaMerged.cpp b/Detector/DetSegmentation/src/FCCSWGridModuleThetaMerged.cpp
index 6c164ddf..8feb2681 100644
--- a/Detector/DetSegmentation/src/FCCSWGridModuleThetaMerged.cpp
+++ b/Detector/DetSegmentation/src/FCCSWGridModuleThetaMerged.cpp
@@ -18,6 +18,7 @@ FCCSWGridModuleThetaMerged::FCCSWGridModuleThetaMerged(const std::string& cellEn
   registerParameter("mergedCells_Theta", "Numbers of merged cells in theta per layer", m_mergedCellsTheta, std::vector<int>());
   registerParameter("mergedModules", "Numbers of merged modules per layer", m_mergedModules, std::vector<int>());
   GetNModulesFromGeom();
+  GetNLayersFromGeom();
 }
 
 FCCSWGridModuleThetaMerged::FCCSWGridModuleThetaMerged(const BitFieldCoder* decoder) : GridTheta(decoder) {
@@ -29,9 +30,9 @@ FCCSWGridModuleThetaMerged::FCCSWGridModuleThetaMerged(const BitFieldCoder* deco
   registerIdentifier("identifier_layer", "Cell ID identifier for layer", m_layerID, "layer");
   registerIdentifier("identifier_module", "Cell ID identifier for module", m_moduleID, "module");
   registerParameter("mergedCells_Theta", "Numbers of merged cells in theta per layer", m_mergedCellsTheta, std::vector<int>());
-  registerParameter("mergedModules", "Numbers of merged cells in phi per layer", m_mergedModules, std::vector<int>());
-  registerParameter("nModules", "Number of modules", m_nModules, 1545);
+  registerParameter("mergedModules", "Numbers of merged modules per layer", m_mergedModules, std::vector<int>());
   GetNModulesFromGeom();
+  GetNLayersFromGeom();
 }
 
 void FCCSWGridModuleThetaMerged::GetNModulesFromGeom() {
@@ -45,6 +46,19 @@ void FCCSWGridModuleThetaMerged::GetNModulesFromGeom() {
   }
   std::cout << "Number of modules read from detector metadata and used in readout class: " << m_nModules << std::endl;
 }
+
+void FCCSWGridModuleThetaMerged::GetNLayersFromGeom() {
+  dd4hep::Detector* dd4hepgeo = &(dd4hep::Detector::getInstance());
+  try {
+    m_nLayers = dd4hepgeo->constant<int>("ECalBarrelNumLayers");
+  }
+  catch(...) {
+    std::cout << "Number of layers not found in detector metadata, exiting..." << std::endl;
+    exit(1);
+  }
+  std::cout << "Number of layers read from detector metadata and used in readout class: " << m_nLayers << std::endl;
+}
+
 /// determine the local position based on the cell ID
 Vector3D FCCSWGridModuleThetaMerged::position(const CellID& cID) const {