[FOCAL-87] Implementing a new design for the HCAL (#14367)

* Implementing a new design for the HCAL.

* Formatting

Author: hahassan7

Detectors/FOCAL/base/include/FOCALBase/Geometry.h

@@ -59,6 +59,12 @@ class VirtualSegment
 class Geometry
 {
  public:
+  enum HCALDesgin {
+    Sandwich = 0,  // Sandwich HCAL design
+    Spaghetti = 1, // Spaghetti design
+    Sheets = 2     // Sheets design
+  };
+
   Geometry() = default;
   Geometry(Geometry* geo);
   Geometry(const Geometry& geo) = default;
@@ -134,6 +140,8 @@ class Geometry
   float getMiddleTowerOffset() const { return mGlobal_Middle_Tower_Offset; }
   bool getInsertFrontPadLayers() const { return mInsertFrontPadLayers; }
   bool getInsertHCalReadoutMaterial() const { return mInsertFrontHCalReadoutMaterial; }
+  float getHCALPitchSize() const { return mGlobal_HCAL_Pitch_Size; }
+  float getHCALBeamPipeHoleSize() const { return mGlobal_HCAL_BeamPipeHole_Size; }
 
   float getDetectorOpeningRight() const { return mGlobal_DetectorOpening_Right; }
   float getDetectorOpeningLeft() const { return mGlobal_DetectorOpening_Left; }
@@ -158,7 +166,7 @@ class Geometry
   void setUpLayerSegmentMap();
   void setUpTowerWaferSize();
 
-  bool getUseHCALSandwich() { return mUseSandwichHCAL; }
+  HCALDesgin getHCALDesign() const { return mHCALDesign; }
 
  protected:
   std::vector<Composition> mGeometryComposition;
@@ -206,9 +214,11 @@ class Geometry
   std::string mGlobal_Gap_Material;         // gap filling material NOTE: currently not used
 
   float mGlobal_HCAL_Tower_Size = 0.0;
-  int mGlobal_HCAL_Tower_NX = 0; // Number of HCAL towers on X
-  int mGlobal_HCAL_Tower_NY = 0; // Number of HCAL towers on Y
-  bool mUseSandwichHCAL = false;
+  int mGlobal_HCAL_Tower_NX = 0;              // Number of HCAL towers on X
+  int mGlobal_HCAL_Tower_NY = 0;              // Number of HCAL towers on Y
+  float mGlobal_HCAL_Pitch_Size = 0.0;        // Distance between two fibers
+  float mGlobal_HCAL_BeamPipeHole_Size = 0.0; // beam pipe hole size in HCAL
+  HCALDesgin mHCALDesign = Sandwich;          // HCAL design type
 
   float mGlobal_FOCAL_Z0 = 0.0;
 

Detectors/FOCAL/base/src/Geometry.cxx

@@ -366,6 +366,11 @@ void Geometry::setParameters(std::string geometryfile)
         LOG(debug) << "The size of the HCAL readout tower will be : " << mGlobal_HCAL_Tower_Size;
       }
 
+      if (command.find("HCAL_PITCH_SIZE") != std::string::npos) {
+        mGlobal_HCAL_Pitch_Size = std::stof(tokens[1]);
+        LOG(debug) << "The distance between fibers is : " << mGlobal_HCAL_Pitch_Size;
+      }
+
       if (command.find("HCAL_TOWER_NX") != std::string::npos) {
         mGlobal_HCAL_Tower_NX = std::stoi(tokens[1]);
         LOG(debug) << "The number of the HCAL readout towers in X will be : " << mGlobal_HCAL_Tower_NX;
@@ -376,6 +381,11 @@ void Geometry::setParameters(std::string geometryfile)
         LOG(debug) << "The number of the HCAL readout towers in Y will be : " << mGlobal_HCAL_Tower_NY;
       }
 
+      if (command.find("HCAL_BEAMPIPE") != std::string::npos) {
+        mGlobal_HCAL_BeamPipeHole_Size = std::stof(tokens[1]);
+        LOG(debug) << "The HCAL beam pipe openning : " << mGlobal_HCAL_BeamPipeHole_Size;
+      }
+
       if (command.find("PIX_OffsetX") != std::string::npos) {
         mGlobal_PIX_OffsetX = std::stof(tokens[1]);
         LOG(debug) << "Pixel offset from the beam pipe will be: " << mGlobal_PIX_OffsetX;
@@ -429,11 +439,6 @@ void Geometry::setParameters(std::string geometryfile)
       if (command.find("NUMBER_OF_HCAL_LAYERS") != std::string::npos) {
         mNHCalLayers = std::stoi(tokens[1]);
         LOG(debug) << "Number of HCAL layers " << mNHCalLayers;
-        if (mNHCalLayers == 1) {
-          mUseSandwichHCAL = false;
-        } else {
-          mUseSandwichHCAL = true;
-        }
       }
 
       if (command.find("NUMBER_OF_SEGMENTS") != std::string::npos) {
@@ -640,6 +645,14 @@ void Geometry::setParameters(std::string geometryfile)
   mHCalLayerThickness = center_z;
   center_z = 0;
 
+  if (mNHCalLayers == 1 && hHCal > 2) {
+    mHCALDesign = Geometry::HCALDesgin::Sheets;
+  } else if (mNHCalLayers == 1 && hHCal == 2) {
+    mHCALDesign = Geometry::HCALDesgin::Spaghetti;
+  } else {
+    mHCALDesign = Geometry::HCALDesgin::Sandwich;
+  }
+
   mFrontMatterLayerThickness = center_z;
   LOG(debug) << " end of SetParameters ";
 }
@@ -702,27 +715,34 @@ std::tuple<double, double, double> Geometry::getGeoTowerCenter(int tower, int se
     int ix = id % nCols;
     int iy = id / nRows;
 
-    if (mUseSandwichHCAL) {
-      float padSize = mVirtualSegmentComposition[segment].mPadSize;
-      double hCALsizeX = nCols * padSize;
-      double hCALsizeY = nRows * padSize;
-      x = ix * padSize + 0.5 * padSize - 0.5 * hCALsizeX;
-      y = iy * padSize + 0.5 * padSize - 0.5 * hCALsizeY;
-    } else {
-      nCols = std::floor(getFOCALSizeX() / getHCALTowerSize() + 0.001) + 1;
-      nRows = std::floor(getFOCALSizeY() / getHCALTowerSize() + 0.001);
-      ix = id % nCols;
-      iy = id / nRows;
-      double beamPipeRadius = 3.6;                                 // in cm   TODO: check if this is OK
-      double towerHalfDiag = std::sqrt(2) * 0.5 * getTowerSizeX(); // tower half diagonal
-      double minRadius = beamPipeRadius + towerHalfDiag;
-
-      float towerSize = getHCALTowerSize() / 7; // To be set from outside (number of channels on x & y)
-      y = iy * towerSize + 0.5 * towerSize - 0.5 * towerSize * nRows;
-      x = ix * towerSize + 0.5 * towerSize - 0.5 * towerSize * nCols;
-      if (y < minRadius && y > -minRadius) {
-        x = int(x) <= 0 ? x - (minRadius - towerSize) : x + (minRadius - towerSize);
+    switch (mHCALDesign) {
+      case HCALDesgin::Sandwich: {
+        float padSize = mVirtualSegmentComposition[segment].mPadSize;
+        double hCALsizeX = nCols * padSize;
+        double hCALsizeY = nRows * padSize;
+
+        x = ix * padSize + 0.5 * padSize - 0.5 * hCALsizeX;
+        y = iy * padSize + 0.5 * padSize - 0.5 * hCALsizeY;
+        break;
       }
+      case HCALDesgin::Spaghetti: {
+        float towerSize = getHCALTowerSize() / 7; // To be set from outside (number of channels on x & y)
+        y = iy * towerSize + 0.5 * towerSize - 0.5 * towerSize * nRows;
+        x = ix * towerSize + 0.5 * towerSize - 0.5 * towerSize * nCols;
+        break;
+      }
+      case HCALDesgin::Sheets: {
+        Composition comp1 = mHCalCompositionBase[0];
+        Composition comp2 = mHCalCompositionBase[2];
+        double hCALsizeX = comp1.sizeX() * 2;                                        // Size of two sheet in X
+        double hCALsizeY = getHCALTowersInY() * (comp1.sizeY() + comp2.sizeY()) * 2; // To be set in a better way
+
+        x = ix * hCALsizeX / getHCALTowersInX() + 0.5 * hCALsizeX / getHCALTowersInX() - 0.5 * hCALsizeX;
+        y = iy * hCALsizeY / getHCALTowersInY() + 0.5 * hCALsizeY / getHCALTowersInY() - 0.5 * hCALsizeY;
+        break;
+      }
+      default:
+        break;
     }
   }
 
@@ -1118,12 +1138,41 @@ std::tuple<bool, int, int, int, int> Geometry::getVirtualInfo(double x, double y
       x = x < 0 ? x - 0.001 : x + 0.001;
       y = y < 0 ? y - 0.001 : y + 0.001;
     }
-    if (!mUseSandwichHCAL) {
-      row = (int)((y + hCALsizeY / 2) / (towerSize / 7));
-      col = (int)((x + hCALsizeX / 2) / (towerSize / 7));
-    } else {
-      row = (int)((y + hCALsizeY / 2) / (towerSize));
-      col = (int)((x + hCALsizeX / 2) / (towerSize));
+
+    switch (mHCALDesign) {
+      case HCALDesgin::Sandwich: {
+        row = (int)((y + hCALsizeY / 2) / (towerSize));
+        col = (int)((x + hCALsizeX / 2) / (towerSize));
+        break;
+      }
+      case HCALDesgin::Spaghetti: {
+        row = (int)((y + hCALsizeY / 2) / (towerSize / 7));
+        col = (int)((x + hCALsizeX / 2) / (towerSize / 7));
+        break;
+      }
+      case HCALDesgin::Sheets: {
+        Composition comp1 = mHCalCompositionBase[0];
+        Composition comp2 = mHCalCompositionBase[2];
+        double hCALsizeX = comp1.sizeX() * 2;                                        // Size of two sheet in X
+        double hCALsizeY = getHCALTowersInY() * (comp1.sizeY() + comp2.sizeY()) * 2; // To be set in a better way
+
+        if (y < getHCALBeamPipeHoleSize() / 2 && y > -getHCALBeamPipeHoleSize() / 2) {
+          if (x < 0) {
+            x += 1.0; // remove the offset around the beam pipe
+          } else {
+            x -= 1.0; // remove the offset around the beam pipe
+          }
+        }
+
+        row = (int)((y + hCALsizeY / 2) / (hCALsizeY / getHCALTowersInY()));
+        if (x > 0) {
+          x = x - 0.15 - 0.06;
+        }
+        col = (int)((x - 0.15 + hCALsizeX / 2) / ((hCALsizeX - 0.15 * 2 - 0.06 * 2) / getHCALTowersInX()));
+        break;
+      }
+      default:
+        break;
     }
   } else {
     row = (int)((y + getFOCALSizeY() / 2) / mVirtualSegmentComposition[segment].mPadSize);
@@ -1150,12 +1199,29 @@ std::tuple<bool, double, double, double> Geometry::getXYZFromColRowSeg(int col,
     double hCALsizeX = getHCALTowersInX() * towerSize;
     double hCALsizeY = getHCALTowersInY() * towerSize;
 
-    if (!mUseSandwichHCAL) {
-      y = -1 * hCALsizeY / 2 + ((float)row + 0.5) * (towerSize / 7);
-      x = -1 * hCALsizeX / 2 + ((float)col + 0.5) * (towerSize / 7);
-    } else {
-      y = -1 * hCALsizeY / 2 + ((float)row + 0.5) * (towerSize);
-      x = -1 * hCALsizeX / 2 + ((float)col + 0.5) * (towerSize);
+    switch (mHCALDesign) {
+      case HCALDesgin::Sandwich: {
+        y = -1 * hCALsizeY / 2 + ((float)row + 0.5) * (towerSize);
+        x = -1 * hCALsizeX / 2 + ((float)col + 0.5) * (towerSize);
+        break;
+      }
+      case HCALDesgin::Spaghetti: {
+        y = -1 * hCALsizeY / 2 + ((float)row + 0.5) * (towerSize / 7);
+        x = -1 * hCALsizeX / 2 + ((float)col + 0.5) * (towerSize / 7);
+        break;
+      }
+      case HCALDesgin::Sheets: {
+        Composition comp1 = mHCalCompositionBase[0];
+        Composition comp2 = mHCalCompositionBase[2];
+        double hCALsizeX = comp1.sizeX() * 2;                                        // Size of two sheet in X
+        double hCALsizeY = getHCALTowersInY() * (comp1.sizeY() + comp2.sizeY()) * 2; // To be set in a better way
+
+        y = -1 * hCALsizeY / 2 + ((float)row + 0.5) * (hCALsizeY / getHCALTowersInY());
+        x = -1 * hCALsizeX / 2 + ((float)col + 0.5) * (hCALsizeX / getHCALTowersInX());
+        break;
+      }
+      default:
+        break;
     }
   } else {
     y = -1 * getFOCALSizeY() / 2 + ((float)row + 0.5) * mVirtualSegmentComposition[segment].mPadSize;
@@ -1191,12 +1257,24 @@ std::tuple<bool, int, int> Geometry::getVirtualNColRow(int segment) const
   nCol = (int)(getFOCALSizeX() / mVirtualSegmentComposition[segment].mPadSize + 0.001);
   nRow = (int)(getFOCALSizeY() / mVirtualSegmentComposition[segment].mPadSize + 0.001);
   if (getVirtualIsHCal(segment)) {
-    if (!mUseSandwichHCAL) {
-      nCol = getHCALTowersInX() * 7; // To be set from outside (number of channels in each tower on x)
-      nRow = getHCALTowersInY() * 7; // To be set from outside (number of channels in each tower on y)
-    } else {
-      nCol = getHCALTowersInX();
-      nRow = getHCALTowersInY();
+    switch (mHCALDesign) {
+      case HCALDesgin::Sandwich: {
+        nCol = getHCALTowersInX();
+        nRow = getHCALTowersInY();
+        break;
+      }
+      case HCALDesgin::Spaghetti: {
+        nCol = getHCALTowersInX() * 7; // To be set from outside (number of channels in each tower on x)
+        nRow = getHCALTowersInY() * 7; // To be set from outside (number of channels in each tower on y)
+        break;
+      }
+      case HCALDesgin::Sheets: {
+        nCol = getHCALTowersInX();
+        nRow = getHCALTowersInY();
+        break;
+      }
+      default:
+        break;
     }
   }
   return {true, nCol, nRow};

Detectors/FOCAL/simulation/geometryFiles/geometry_Sheets.txt

@@ -0,0 +1,86 @@
+# Volume Name            Material     C-X   C-Y    X      Y       Z
+  COMPOSITION_PAD_S0       Alloy       0     0     9.     8.      0.35
+  COMPOSITION_PAD_S1       G10         0     0     9.     8.      0.08
+  COMPOSITION_PAD_S2       SiPad       0     0     9.     8.      0.03
+  COMPOSITION_PAD_S3       G10         0     0     9.     8.      0.08
+  COMPOSITION_PAD_S4       Cu          0     0     9.     8.      0.014
+  COMPOSITION_PAD_S5       Air         0     0     9.     8.      0.296
+  # Replica of above pad layers to 50 layers
+  COMMAND_NUMBER_OF_PAD_LAYERS 18
+# HCAL layers
+# Volume Name              Material     C-X   C-Y     X       Y        Z
+  COMPOSITION_HCAL_S0       CuHCAL       0     0     49.81   0.20     110
+  COMPOSITION_HCAL_S1       Scint        0     0     0.1     0.1      110
+  COMPOSITION_HCAL_S2       CuHCAL       0     0     49.81   0.15     110
+  COMMAND_NUMBER_OF_HCAL_LAYERS 1
+  COMMAND_NUMBER_OF_SEGMENTS 21
+# Strip sectors
+# Volume Name             Material    C-X   C-Y     X       Y         Z
+  COMPOSITION_PIX_S0       Alloy       0     0     3.0     2.74      0.35
+  COMPOSITION_PIX_S1       G10         0     0     3.0     2.74      0.1
+  COMPOSITION_PIX_S2       Si          0     0     3.0     2.74      0.047
+  COMPOSITION_PIX_S3       SiPix       0     0     3.0     2.74      0.003
+  COMPOSITION_PIX_S4       G10         0     0     3.0     2.74      0.1
+  COMPOSITION_PIX_S5       Cu          0     0     3.0     2.74      0.001
+  COMPOSITION_PIX_S6       Air         0     0     3.0     2.74      0.249
+COMMAND_INSERT_PIX_AT_L4
+COMMAND_INSERT_PIX_AT_L9
+# COMMAND_INSERT_STR_AT_L4
+# COMMAND_INSERT_STR_AT_L6
+# COMMAND_INSERT_STR_AT_L8
+#Front Matter definition
+  COMPOSITION_FM_S0        G10         0     0     5.    5.      0.01
+  COMPOSITION_FM_S1        SiStripX    0     0     5.    5.      0.05
+  COMPOSITION_FM_S2        G10         0     0     5.    5.      0.05
+  COMPOSITION_FM_S3        Air         0     0     5.    5.      0.01
+  COMPOSITION_FM_S4        G10         0     0     5.    5.      0.01
+  COMPOSITION_FM_S5        SiStripY    0     0     5.    5.      0.05
+  COMPOSITION_FM_S6        G10         0     0     5.    5.      0.05
+# PIXEL readout
+  GLOBAL_PIX_NX                  15
+  GLOBAL_PIX_NY                  3
+  GLOBAL_PIX_OffsetX             1.2
+  GLOBAL_PIX_OffsetY             0.09
+  GLOBAL_PIX_SKIN                0.004
+#  pxel size in cm
+  COMMAND_PIXEL_READOUT_ON       0.005
+# Pad information
+  GLOBAL_PAD_SIZE_X_Y            1
+  GLOBAL_PAD_NX_NY               8
+  GLOBAL_PAD_NX                  9
+  GLOBAL_PAD_NY                  8
+  GLOBAL_PAD_PPTOL               0.
+  GLOBAL_PAD_SKIN                0.2
+# Global information (TOL:1cm of T-T space filled with "AIR")
+  GLOBAL_PAD_SUPERMODULE_X 5
+  GLOBAL_PAD_SUPERMODULE_Y 1
+  GLOBAL_SUPERMODULE_TOLX 0. Air
+  GLOBAL_SUPERMODULE_TOLY 0. Air
+  GLOBAL_TOWER_TOL        0.    Air
+  GLOBAL_TOWER_TOLX       0.02  Air
+  GLOBAL_TOWER_TOLY       0.8   Al
+  GLOBAL_FOCAL_Z          764.47
+  GLOBAL_Tower_NX      2
+  GLOBAL_Tower_NY      11
+  GLOBAL_MIDDLE_TOWER_OFFSET 5
+  GLOBAL_NSTRIPS	        128
+  GLOBAL_STRIPSIZE_LONG   9.0
+  GLOBAL_STRIPSIZE_WIDTH  0.07
+  GLOBAL_HCAL_PITCH_SIZE  0.4
+  GLOBAL_HCAL_TOWER_NY  72
+  GLOBAL_HCAL_TOWER_NX  62
+  GLOBAL_HCAL_BEAMPIPE  8.4
+# COMMAND_INSERT_FRONT_PAD_LAYERS
+# COMMAND_INSERT_HCAL_READOUT
+# New VIRTUAL settings
+  VIRTUAL_N_SEGMENTS      7
+  #                            N   Start End   PadSize  RelThickness  IsPixel    PixelTreshold [eV]
+  #                                Layer Layer
+  VIRTUAL_SEGMENT_LAYOUT_N0       0    3           1.0           1.0       0                  300000
+  VIRTUAL_SEGMENT_LAYOUT_N1       4    4           0.05          1.0       1                  4000
+  VIRTUAL_SEGMENT_LAYOUT_N2       5    8           1.0           1.0       0                  300000
+  VIRTUAL_SEGMENT_LAYOUT_N3       9    9           0.05          1.0       1                  4000
+  VIRTUAL_SEGMENT_LAYOUT_N4      10   14           1.0           1.0       0                  375000
+  VIRTUAL_SEGMENT_LAYOUT_N5      15   19           1.0           1.0       0                  375000
+  VIRTUAL_SEGMENT_LAYOUT_N6      20   20           0.4           1.0       2                  5000
+#EOF

Detectors/FOCAL/simulation/include/FOCALSimulation/Detector.h

@@ -17,6 +17,7 @@
 #include "DetectorsBase/Detector.h"
 #include "FOCALBase/Hit.h"
 #include "FOCALBase/Geometry.h"
+#include "TGeoManager.h"
 
 class FairVolume;
 
@@ -157,6 +158,17 @@ class Detector : public o2::base::DetImpl<Detector>
 
   virtual void CreateHCALSpaghetti();
   virtual void CreateHCALSandwich();
+  virtual void CreateHCALSheets();
+
+  TGeoVolumeAssembly* CreatePitchAssembly(double Lx = 498.1,
+                                          double Ly1 = 2.0,
+                                          double Ly2 = 1.5,
+                                          double Lz = 1100.0,
+                                          double hole_diameter = 1.1,
+                                          double hole_spacing = 4.0,
+                                          int nholes = 124,
+                                          double fiber_radius = 0.5,
+                                          std::string suffix = "");
 
   /// \brief Generate ECAL geometry
   void CreateECALGeometry();