Support arbitrary cross-section shapes in beamline benchmark overlays

.gitignore

@@ -51,3 +51,4 @@ calorimeters/test/
 # org2py output
 benchmarks/backgrounds/ecal_backwards.py
 benchmarks/ecal_gaps/ecal_gaps.py
+_codeql_detected_source_root

benchmarks/beamline/acceptanceAnalysis.C

@@ -111,6 +111,13 @@ int acceptanceAnalysis( TString inFile             = "/home/simong/EIC/detector_
                     radii.push_back(param.radius);
                 }
                 return radii;
+                }, {"pipeParameters"})
+                .Define("isConeSegment",[](const ROOT::VecOps::RVec<volParams>& params) {
+                ROOT::VecOps::RVec<bool> cones;
+                for (const auto& param : params) {
+                    cones.push_back(param.isConeSegment);
+                }
+                return cones;
                 }, {"pipeParameters"});
 
 
@@ -139,6 +146,8 @@ int acceptanceAnalysis( TString inFile             = "/home/simong/EIC/detector_
 
 
     std::map<TString,ROOT::RDF::RResultPtr<double>> pipeRadii;
+    std::map<TString,ROOT::RDF::RResultPtr<bool>> pipeIsConeSegment;
+    std::map<TString,ROOT::RDF::RResultPtr<volParams>> pipeParams;
     std::map<TString,double> filterEntries;
     std::map<TString,double> filterAcceptanceIntegral;
 
@@ -150,7 +159,9 @@ int acceptanceAnalysis( TString inFile             = "/home/simong/EIC/detector_
         name += str_i;
         auto filterDF = d1.Define("xposf","xpos[pipeID=="+str_i+"]")
                           .Define("yposf","ypos[pipeID=="+str_i+"]")
-                          .Define("pipeRadiusf","pipeRadius[pipeID=="+str_i+"]");
+                          .Define("pipeRadiusf","pipeRadius[pipeID=="+str_i+"]")
+                          .Define("isConeSegmentf","isConeSegment[pipeID=="+str_i+"]")
+                          .Define("pipeParametersf","pipeParameters[pipeID=="+str_i+"]");
 
 
         TString beamspotName = "Beamspot ID"+str_i+";x offset [cm]; y offset [cm]";
@@ -165,7 +176,9 @@ int acceptanceAnalysis( TString inFile             = "/home/simong/EIC/detector_
         hHistsETheta[name] = extraFilterDF.Histo2D({EThetaName,EThetaName,eBins,4,18,thetaBins,0,0.011},"energy","theta");
 
         //Parameters of the pipe
-        pipeRadii[name]    = filterDF.Max("pipeRadiusf");        
+        pipeRadii[name]    = filterDF.Max("pipeRadiusf");
+        pipeIsConeSegment[name] = filterDF.Max("isConeSegmentf");
+        pipeParams[name] = filterDF.Max("pipeParametersf");
         // std::cout << "Pipe ID: " << name << " Radius: " << pipeRadii[name] << " " << filterDF.Min("pipeRadiusf").GetValue() << std::endl;
 
     }
@@ -175,17 +188,25 @@ int acceptanceAnalysis( TString inFile             = "/home/simong/EIC/detector_
     cXY->Divide(4,2);
     int i=1;
     for(auto [name,h] : hHistsxy){
-        // Get the pipe radius for this histogram
+        // Get the pipe parameters for this histogram
         auto pipeRadius = pipeRadii[name].GetValue();
+        auto isCone = pipeIsConeSegment[name].GetValue();
+        auto params = pipeParams[name].GetValue();
         cXY->cd(i++);
 
         h->Draw("col");
-        //Draw cicle
-        TEllipse *circle = new TEllipse(0,0,pipeRadius);
-        circle->SetLineColor(kRed);
-        circle->SetLineWidth(2);
-        circle->SetFillStyle(0);
-        circle->Draw("same");
+
+        // Only draw circle overlay if the shape is a cone segment
+        if (isCone && pipeRadius > 0) {
+            TEllipse *circle = new TEllipse(0,0,pipeRadius);
+            circle->SetLineColor(kRed);
+            circle->SetLineWidth(2);
+            circle->SetFillStyle(0);
+            circle->Draw("same");
+        } else if (!params.shapeOutline.empty()) {
+            // Draw arbitrary shape outline for non-cone segments
+            drawShapeOutline(params.shapeOutline, kRed, 2);
+        }
 
     }
 

benchmarks/beamline/beamlineAnalysis.C

@@ -178,6 +178,7 @@ int beamlineAnalysis(   TString inFile          = "/home/simong/EIC/detector_ben
     std::map<TString,double> pipeZPos;
     std::map<TString,double> pipeRotation;
     std::map<TString,bool> pipeIsConeSegment;
+    std::map<TString,volParams> pipeParams;
 
     // Queue up all actions
     auto xmin_ptr  = d1.Min("xpos");
@@ -210,6 +211,7 @@ int beamlineAnalysis(   TString inFile          = "/home/simong/EIC/detector_ben
                           .Define("ymomf","ymom[pipeID=="+std::to_string(i)+"]")
                           .Define("pipeRadiusf","pipeRadius[pipeID=="+std::to_string(i)+"]")
                           .Define("isConeSegmentf","isConeSegment[pipeID=="+std::to_string(i)+"]")
+                          .Define("pipeParametersf","pipeParameters[pipeID=="+std::to_string(i)+"]")
                           .Define("xdetf","xdet[pipeID=="+std::to_string(i)+"]")
                           .Define("zdetf","zdet[pipeID=="+std::to_string(i)+"]")
                           .Define("rotationf","rotation[pipeID=="+std::to_string(i)+"]");
@@ -283,6 +285,7 @@ int beamlineAnalysis(   TString inFile          = "/home/simong/EIC/detector_ben
         pipeZPos[name]     = filterDF.Max("zdetf").GetValue();
         pipeRotation[name] = filterDF.Max("rotationf").GetValue();
         pipeIsConeSegment[name] = filterDF.Max("isConeSegmentf").GetValue();
+        pipeParams[name] = filterDF.Max("pipeParametersf").GetValue();
 
         //Fit gaussian to the x, y, px and py distributions
         auto xhist = hHistsxy[name]->ProjectionX();
@@ -343,6 +346,9 @@ int beamlineAnalysis(   TString inFile          = "/home/simong/EIC/detector_ben
             circle->SetLineWidth(2);
             circle->SetFillStyle(0);
             circle->Draw("same");
+        } else if (!pipeParams[name].shapeOutline.empty()) {
+            // Draw arbitrary shape outline for non-cone segments
+            drawShapeOutline(pipeParams[name].shapeOutline, kRed, 2);
         }
 
         // Add zoomed version in the top-right corner

benchmarks/beamline/shared_functions.h

@@ -6,6 +6,9 @@
 #include "DD4hep/VolumeManager.h"
 #include "DD4hep/DetElement.h"
 #include "TFile.h"
+#include "TGeoShape.h"
+#include "TGeoBBox.h"
+#include "TPolyLine.h"
 
 using RVecHits = ROOT::VecOps::RVec<edm4hep::SimTrackerHitData>;
 using namespace dd4hep;
@@ -104,13 +107,92 @@ struct globalToLocal{
   dd4hep::VolumeManager volumeManager;
 };
 
+//-----------------------------------------------------------------------------------------
+// Helper function to extract shape outline points for arbitrary TGeo shapes
+//-----------------------------------------------------------------------------------------
+std::vector<std::pair<double, double>> extractShapeOutline(dd4hep::Solid solid, double z = 0.0, int nPoints = 100) {
+  std::vector<std::pair<double, double>> outline;
+
+  if (!solid.isValid()) {
+    return outline;
+  }
+
+  TGeoShape* geoShape = solid.ptr();
+  std::string shapeName = geoShape->IsA()->GetName();
+
+  // For ConeSegment, return circular outline
+  if (shapeName == "TGeoConeSeg") {
+    dd4hep::ConeSegment cone = solid;
+    double rmax = cone.rMax1();
+    for (int i = 0; i <= nPoints; i++) {
+      double angle = 2.0 * M_PI * i / nPoints;
+      outline.push_back({rmax * std::cos(angle), rmax * std::sin(angle)});
+    }
+    return outline;
+  }
+
+  // For other shapes, sample points on the boundary
+  // This is a general approach that works for intersection solids and other complex shapes
+  TGeoBBox* bbox = dynamic_cast<TGeoBBox*>(geoShape);
+  if (bbox) {
+    // Get bounding box dimensions
+    double dx = bbox->GetDX();
+    double dy = bbox->GetDY();
+    // For a box, trace the outline
+    outline.push_back({-dx, -dy});
+    outline.push_back({ dx, -dy});
+    outline.push_back({ dx,  dy});
+    outline.push_back({-dx,  dy});
+    outline.push_back({-dx, -dy});
+    return outline;
+  }
+
+  // For composite shapes (like intersection solids), sample the boundary
+  // by testing points in a grid and finding those near the surface
+  double maxR = 10.0; // Maximum radius to check (in cm)
+  double step = 0.1;  // Step size for sampling
+
+  // Sample points in polar coordinates to find the shape boundary
+  for (int i = 0; i <= nPoints; i++) {
+    double angle = 2.0 * M_PI * i / nPoints;
+    double cosAngle = std::cos(angle);
+    double sinAngle = std::sin(angle);
+
+    // Binary search to find the boundary at this angle
+    double rMin = 0.0;
+    double rMax = maxR;
+    double r = rMax / 2.0;
+
+    for (int iter = 0; iter < 20; iter++) { // 20 iterations for binary search
+      double x = r * cosAngle;
+      double y = r * sinAngle;
+      double point[3] = {x, y, z};
+
+      if (geoShape->Contains(point)) {
+        rMin = r;
+      } else {
+        rMax = r;
+      }
+      r = (rMin + rMax) / 2.0;
+    }
+
+    // If we found a valid boundary point
+    if (rMin > 0.01) { // Threshold to avoid noise
+      outline.push_back({rMin * cosAngle, rMin * sinAngle});
+    }
+  }
+
+  return outline;
+}
+
 struct volParams{
   double radius;
   double xPos;
   double yPos;
   double zPos;
   double rotation;
   bool isConeSegment;
+  std::vector<std::pair<double, double>> shapeOutline; // (x, y) coordinates in local frame
 };
 
 // Functor to get the volume element parameters from the CellID
@@ -137,13 +219,17 @@ struct getVolumeParametersFromCellID {
         dd4hep::ConeSegment cone = solid;
         radius = cone.rMax1();
       }
+      // Extract shape outline for arbitrary shapes
+      std::vector<std::pair<double, double>> outline = extractShapeOutline(solid);
+
       volParams params{
         radius,
         translation[0],
         translation[1],
         translation[2],
         rotationAngleY,
-        isCone
+        isCone,
+        outline
       };
       result.push_back(params);
     }
@@ -155,6 +241,32 @@ struct getVolumeParametersFromCellID {
   dd4hep::VolumeManager volumeManager;
 };
 
+//-----------------------------------------------------------------------------------------
+// Helper function to draw shape outline on a canvas
+//-----------------------------------------------------------------------------------------
+void drawShapeOutline(const std::vector<std::pair<double, double>>& outline, 
+                      int lineColor = kRed, 
+                      int lineWidth = 2) {
+  if (outline.empty()) {
+    return;
+  }
+
+  int nPoints = outline.size();
+  std::vector<double> x(nPoints);
+  std::vector<double> y(nPoints);
+
+  for (int i = 0; i < nPoints; i++) {
+    x[i] = outline[i].first;
+    y[i] = outline[i].second;
+  }
+
+  TPolyLine* poly = new TPolyLine(nPoints, x.data(), y.data());
+  poly->SetLineColor(lineColor);
+  poly->SetLineWidth(lineWidth);
+  poly->SetFillStyle(0);
+  poly->Draw("same");
+}
+
 TH1F* CreateFittedHistogram(const std::string& histName, 
   const std::string& histTitle, 
   const std::map<TString, double>& valueMap,