Fiber DRC segmentation and SD action update for downstream digitization and reconstruction

FCCee/IDEA/compact/IDEA_o1_v03/FiberDualReadoutCalo_o1_v01.xml

@@ -594,6 +594,12 @@
       </sipmDim>
       <structure>
         <dim distance="1.5*mm" dx="1.0*mm"/>
+        <!--parameters for neighborhood finding algorithm-->
+        <!--rmax1/rmax2: radius of the neighborhood in units of fiber-fiber distance in barrel/endcap-->
+        <!--x1/x2: margin to incorporate satelite cells at the edge of the tower in barrel/endcap-->
+        <!--delta: flag to remove different type of cells in the neighborhood-->
+        <!--i.e. if set to 0 then Cherenkov and scint fibers will be grouped together-->
+        <neighborhood rmax1="1.5" rmax2="1.5" x1="2" x2="2" delta="0"/>
         <cladC name="cladC" rmax="0.5*mm" rmin="0.49*mm" material="FluorinatedPolymer" vis="DRCCladVis"/>
         <coreC name="coreC" rmax="0.5*mm" rmin="0.49*mm" material="PMMA" vis="DRCCerenVis"/> <!--also cladS-->
         <coreS name="coreS" rmax="0.5*mm" rmin="0.485*mm" material="DR_Polystyrene" vis="DRCScintVis"/>

FCCee/IDEA/compact/README.md

@@ -30,6 +30,8 @@ Febraury 2025: Added surface plugins for both muon-system and preshower, with ad
 
 March 2025: The SiPM & optical filter part of the fiber DRC SD action are now moved to the k4RecCalorimeter (now part of the SiPM simulation)
 
+July 2025: Added the fiber DRC neighborhood finding algorithm with configure size definition for topological clustering in downstream.
+
 IDEA_o2_v01
 ------------
 

detector/calorimeter/README.md

@@ -62,6 +62,7 @@ This sub-detector makes full 4-pi monolithic fiber dual-readout calorimeter.
 Inside the single tower (trapezoidal copper absorber), two types of optical fibers (Cherenkov and scintillation) are implemented. The readout (SiPM) is attached at the rear side of the tower. The tower is repeated in both eta and phi direction to cover both barrel and endcap region.
 
 Added an extension to (`LayeredCalorimeterData`) to store the barrel and endcap rmin, rmax, zmin, zmax.
+The definition of the topological neighborhood of cells is configurable via the xml file (function itself is defined in the `GridDRcalo_k4geo`).
 
 ## Capillary tube dual-readout (subdirectory `dual-readout-tubes`)
 

detector/calorimeter/dual-readout/src/DRconstructor.cpp

@@ -18,7 +18,6 @@ ddDRcalo::DRconstructor::DRconstructor(xml_det_t& x_det)
   fVis = false;
   fNumx = 0;
   fNumy = 0;
-  fFiberCoords.reserve(100000);
   fFiberEnvVec.reserve(4000);
   fFiberCoreCVec.reserve(4000);
   fFiberCoreSVec.reserve(4000);
@@ -52,6 +51,7 @@ void ddDRcalo::DRconstructor::construct() {
   // since reflecting intersection/subtraction solid is currently not possible
   // this causes a discontinuity of the neighboring cells at eta=0
   // we handle this explicitly in the segmentation
+  // be aware that it needs to be synchronized with GridDRcalo_k4geo::position()
   if (fX_det.reflect()) {
     auto refl_pos = dd4hep::Transform3D(dd4hep::RotationX(M_PI), dd4hep::Position(0, 0, -(fX_worldTube.height() / 2.)));
     dd4hep::PlacedVolume PlacedAssemblyTubeVol_refl = fExperimentalHall->placeVolume(AssemblyTubeVol, 1, refl_pos);
@@ -109,9 +109,6 @@ void ddDRcalo::DRconstructor::implementTowers(xml_comp_t& x_theta, dd4hep::DDSeg
       sipmWaferVol.setSensitiveDetector(*fSensDet);
     }
 
-    // Remove sipmLayer, clear fFiberCoords instead of implementSipms()
-    fFiberCoords.clear();
-
     for (int nPhi = 0; nPhi < x_theta.nphi(); nPhi++) {
       placeAssembly(param, AssemblyBoxVol, towerVol, sipmWaferVol, towerNo, nPhi);
     }
@@ -189,6 +186,9 @@ void ddDRcalo::DRconstructor::implementFibers(xml_comp_t& x_theta, dd4hep::Volum
   double norm1[3] = {0., 0., 0.}, norm2[3] = {0., 0., 0.}, norm3[3] = {0., 0., 0.}, norm4[3] = {0., 0., 0.};
   getNormals(rootTrap, numxBl2, z1, norm1, norm2, norm3, norm4);
 
+  // std::map containing the length of short fibers
+  dd4hep::DDSegmentation::DRparamBase_k4geo::shortFibers shortFibers(towerHeight);
+
   for (int row = 0; row < fNumy; row++) {
     for (int column = 0; column < fNumx; column++) {
       auto localPosition = fSegmentation->localPosition(fNumx, fNumy, column, row);
@@ -201,7 +201,6 @@ void ddDRcalo::DRconstructor::implementFibers(xml_comp_t& x_theta, dd4hep::Volum
 
           if (check) {
             implementFiber(fullBoxVol, pos, column, row);
-            fFiberCoords.push_back(std::make_pair(column, row));
           }
         }
       } else {
@@ -241,14 +240,15 @@ void ddDRcalo::DRconstructor::implementFibers(xml_comp_t& x_theta, dd4hep::Volum
         if (checkContained(rootTrap, pos, towerHeight / 2. - fiberLen)) {
           dd4hep::Position centerPos(pos.x(), pos.y(), centerZ);
           implementFiber(towerVol, centerPos, column, row, fiberLen);
-          fFiberCoords.push_back(std::make_pair(column, row));
+          shortFibers.addShortFibers(row, column, fiberLen);
         }
-      }
-    }
-  }
+      } // outside tower
+    } // col
+  } // row
 
   // store rows & columns of full length fibers to the segmentation
   param->SetFullLengthFibers(rmin, rmax, cmin, cmax);
+  param->SetShortFibers(shortFibers);
 }
 
 // Remove cap (mirror or black paint in front of the fiber)

detector/calorimeter/dual-readout/src/FiberDualReadoutCalo_o1_v01.cpp

@@ -31,6 +31,7 @@ static dd4hep::Ref_t create_detector(dd4hep::Detector& description, xml_h xmlEle
   xml_comp_t x_dim(x_structure.child(_Unicode(dim)));
   xml_comp_t x_sipmDim(x_det.child(_Unicode(sipmDim)));
   xml_comp_t x_worldTube(x_structure.child(_Unicode(worldTube)));
+  xml_comp_t x_neighbor(x_structure.child(_Unicode(neighborhood)));
 
   dd4hep::OpticalSurfaceManager surfMgr = description.surfaceManager();
   dd4hep::OpticalSurface sipmSurfProp = surfMgr.opticalSurface("/world/" + name + "#SiPMSurf");
@@ -40,18 +41,23 @@ static dd4hep::Ref_t create_detector(dd4hep::Detector& description, xml_h xmlEle
       dynamic_cast<dd4hep::DDSegmentation::GridDRcalo_k4geo*>(sensDet.readout().segmentation().segmentation());
   segmentation->setGridSize(x_dim.distance());
   segmentation->setSipmSize(x_dim.dx());
+  segmentation->setRemoveDifferentCh(static_cast<bool>(x_neighbor.delta()));
 
   auto paramBarrel = segmentation->paramBarrel();
   paramBarrel->SetInnerX(x_barrel.rmin());
   paramBarrel->SetTowerH(x_barrel.height());
   paramBarrel->SetNumZRot(x_barrel.nphi());
   paramBarrel->SetSipmHeight(x_sipmDim.height());
+  paramBarrel->SetNeighborSize(x_neighbor.rmax1());
+  paramBarrel->SetMargin(x_neighbor.x1());
 
   auto paramEndcap = segmentation->paramEndcap();
   paramEndcap->SetInnerX(x_endcap.rmin());
   paramEndcap->SetTowerH(x_endcap.height());
   paramEndcap->SetNumZRot(x_endcap.nphi());
   paramEndcap->SetSipmHeight(x_sipmDim.height());
+  paramEndcap->SetNeighborSize(x_neighbor.rmax2());
+  paramEndcap->SetMargin(x_neighbor.x2());
 
   auto constructor = DRconstructor(x_det);
   constructor.setExpHall(&experimentalHall);

detectorSegmentations/include/detectorSegmentations/DRparamBase_k4geo.h

@@ -50,8 +50,8 @@ namespace DDSegmentation {
     dd4hep::Transform3D GetAssembleTransform3D(int numPhi);
     dd4hep::Transform3D GetSipmTransform3D(int numPhi);
 
-    int signedTowerNo(int unsignedTowerNo) { return fIsRHS ? unsignedTowerNo : -unsignedTowerNo - 1; }
-    int unsignedTowerNo(int signedTowerNo) { return signedTowerNo >= 0 ? signedTowerNo : -signedTowerNo - 1; }
+    int signedTowerNo(int unsignedTowerNo) const { return fIsRHS ? unsignedTowerNo : -unsignedTowerNo - 1; }
+    int unsignedTowerNo(int signedTowerNo) const { return signedTowerNo >= 0 ? signedTowerNo : -signedTowerNo - 1; }
 
     virtual void SetDeltaThetaByTowerNo(int, int) {}
     virtual void SetThetaOfCenterByTowerNo(int, int) {}
@@ -63,6 +63,13 @@ namespace DDSegmentation {
     int GetCurrentTowerNum() { return fCurrentTowerNum; }
     void SetCurrentTowerNum(int numEta) { fCurrentTowerNum = numEta; }
 
+    // the size of the neighborhood (radius) in units of fiber-fiber distance (i.e. grid size)
+    void SetNeighborSize(double nsize) { fNeighborSize = nsize; }
+    double GetNeighborSize() const { return fNeighborSize; }
+    // the margin in units of fiber-fiber distance when incorporating satelite cells at the edge of the tower
+    void SetMargin(int margin) { fMargin = margin; }
+    int GetMargin() const { return fMargin; }
+
     virtual void init() {};
     void filled() { fFilled = true; }
     void finalized() { fFinalized = true; }
@@ -84,9 +91,30 @@ namespace DDSegmentation {
       int cmax; // max n_column
     };
 
-    fullLengthFibers GetFullLengthFibers(int numEta) { return fFullLengthFibers.at(unsignedTowerNo(numEta)); }
+    fullLengthFibers GetFullLengthFibers(int numEta) const { return fFullLengthFibers.at(unsignedTowerNo(numEta)); }
     void SetFullLengthFibers(int rmin, int rmax, int cmin, int cmax);
 
+    // length of fibers that don't have full length
+    // it is a map containing (row, col) as a key and the length of the fiber as a value
+    // and each tower in eta has one of it
+    struct shortFibers {
+    public:
+      shortFibers(const double towerH) : fTowerH(towerH) {}
+      ~shortFibers() = default;
+
+      void addShortFibers(const int row, const int col, const double len) {
+        m_fiberLengths_.insert(std::make_pair(std::make_pair(row, col), len));
+      }
+      double retrieveFiberLength(const int row, const int col) const;
+
+    private:
+      std::map<std::pair<int, int>, double> m_fiberLengths_;
+      const double fTowerH;
+    };
+
+    shortFibers GetShortFibers(int numEta) const { return fShortFibers.at(unsignedTowerNo(numEta)); }
+    void SetShortFibers(const shortFibers& input);
+
   protected:
     bool fIsRHS;
     double fPhiZRot;
@@ -109,11 +137,16 @@ namespace DDSegmentation {
     double fCurrentOuterHalf;
     double fCurrentOuterHalfSipm;
 
+    // parameters for neighbour finding algorithm
+    double fNeighborSize;
+    int fMargin;
+
     int fTotNum;
     int fCurrentTowerNum;
     std::vector<double> fDeltaThetaVec;
     std::vector<double> fThetaOfCenterVec;
     std::map<int, fullLengthFibers> fFullLengthFibers;
+    std::map<int, shortFibers> fShortFibers;
     bool fFilled;
     bool fFinalized;
   };

detectorSegmentations/include/detectorSegmentations/GridDRcaloHandle_k4geo.h

@@ -34,6 +34,7 @@ class GridDRcalo_k4geo : public GridDRcaloHandle {
   bool operator==(const GridDRcalo_k4geo& seg) const { return m_element == seg.m_element; }
   /// determine the position based on the cell ID
   inline Position position(const CellID& id) const { return Position(access()->implementation->position(id)); }
+  inline Position sipmPosition(const CellID& id) const { return Position(access()->implementation->sipmPosition(id)); }
   inline Position localPosition(const CellID& id) const {
     return Position(access()->implementation->localPosition(id));
   }

detectorSegmentations/include/detectorSegmentations/GridDRcalo_k4geo.h

@@ -17,11 +17,15 @@ namespace DDSegmentation {
     /// destructor
     virtual ~GridDRcalo_k4geo() override;
 
-    //  Determine the global(local) position based on the cell ID.
+    // Determine the global (local) position based on the cell ID.
+    // position of the front end of the fiber
     virtual Vector3D position(const CellID& aCellID) const override;
     Vector3D localPosition(const CellID& aCellID) const;
     Vector3D localPosition(int numx, int numy, int x_, int y_) const;
 
+    // position of the rear end (readout location) of the fiber
+    Vector3D sipmPosition(const CellID& aCellID) const;
+
     virtual CellID cellID(const Vector3D& aLocalPosition, const Vector3D& aGlobalPosition,
                           const VolumeID& aVolumeID) const override;
 
@@ -33,6 +37,9 @@ namespace DDSegmentation {
     void setGridSize(double grid) { fGridSize = grid; }
     void setSipmSize(double sipm) { fSipmSize = sipm; }
 
+    // remove different type of channels in the neighborhood set if set to true
+    void setRemoveDifferentCh(bool rmDiffCh) { fRemoveDifferentCh = rmDiffCh; }
+
     // Get the identifier number of a mother tower in eta or phi direction
     int numEta(const CellID& aCellID) const;
     int numPhi(const CellID& aCellID) const;
@@ -105,7 +112,9 @@ namespace DDSegmentation {
     double fGridSize;
     double fSipmSize;
 
-  private:
+    // switch to remove different type cells from the neighborhood
+    bool fRemoveDifferentCh;
+
     DRparamBarrel_k4geo* fParamBarrel;
     DRparamEndcap_k4geo* fParamEndcap;
   };

detectorSegmentations/src/DRparamBase_k4geo.cpp

@@ -108,5 +108,20 @@ namespace DDSegmentation {
     fFullLengthFibers.insert(std::make_pair(fCurrentTowerNum, fullLengthFibers(rmin, rmax, cmin, cmax)));
   }
 
+  double DRparamBase_k4geo::shortFibers::retrieveFiberLength(const int row, const int col) const {
+    if (m_fiberLengths_.find(std::make_pair(row, col)) == m_fiberLengths_.end())
+      return fTowerH;
+
+    return m_fiberLengths_.at(std::make_pair(row, col));
+  }
+
+  void DRparamBase_k4geo::SetShortFibers(const shortFibers& input) {
+    if (fFilled)
+      throw std::runtime_error(
+          "DRparamBase_k4geo: An attempt to modify the geometry outside the detector construction is forbidden!");
+
+    fShortFibers.insert(std::make_pair(fCurrentTowerNum, input));
+  }
+
 } // namespace DDSegmentation
 } // namespace dd4hep