diff --git a/.gitignore b/.gitignore index 66337cd9..30756b19 100644 --- a/.gitignore +++ b/.gitignore @@ -51,3 +51,4 @@ calorimeters/test/ # org2py output benchmarks/backgrounds/ecal_backwards.py benchmarks/ecal_gaps/ecal_gaps.py +_codeql_detected_source_root diff --git a/benchmarks/beamline/acceptanceAnalysis.C b/benchmarks/beamline/acceptanceAnalysis.C index 800ba024..f11a43dd 100644 --- a/benchmarks/beamline/acceptanceAnalysis.C +++ b/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& params) { + ROOT::VecOps::RVec 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> pipeRadii; + std::map> pipeIsConeSegment; + std::map> pipeParams; std::map filterEntries; std::map 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); + } } diff --git a/benchmarks/beamline/beamlineAnalysis.C b/benchmarks/beamline/beamlineAnalysis.C index 08ae6e33..2fef465c 100644 --- a/benchmarks/beamline/beamlineAnalysis.C +++ b/benchmarks/beamline/beamlineAnalysis.C @@ -178,6 +178,7 @@ int beamlineAnalysis( TString inFile = "/home/simong/EIC/detector_ben std::map pipeZPos; std::map pipeRotation; std::map pipeIsConeSegment; + std::map 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 diff --git a/benchmarks/beamline/shared_functions.h b/benchmarks/beamline/shared_functions.h index b0bf1e8e..b9d3f6b5 100644 --- a/benchmarks/beamline/shared_functions.h +++ b/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; using namespace dd4hep; @@ -104,6 +107,84 @@ struct globalToLocal{ dd4hep::VolumeManager volumeManager; }; +//----------------------------------------------------------------------------------------- +// Helper function to extract shape outline points for arbitrary TGeo shapes +//----------------------------------------------------------------------------------------- +std::vector> extractShapeOutline(dd4hep::Solid solid, double z = 0.0, int nPoints = 100) { + std::vector> 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(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; @@ -111,6 +192,7 @@ struct volParams{ double zPos; double rotation; bool isConeSegment; + std::vector> 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> 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>& outline, + int lineColor = kRed, + int lineWidth = 2) { + if (outline.empty()) { + return; + } + + int nPoints = outline.size(); + std::vector x(nPoints); + std::vector 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& valueMap,