Beamline - Allow detector volumes which are not cone segments for tracking (#221)

Author: simonge

benchmarks/beamline/acceptanceAnalysis.C

@@ -19,7 +19,7 @@
 using RVecS       = ROOT::VecOps::RVec<string>;
 using RNode       = ROOT::RDF::RNode;
 
-int acceptanceAnalysis( TString inFile             = "/scratch/EIC/G4out/beamline/acceptanceTest.edm4hep.root",
+int acceptanceAnalysis( TString inFile             = "/home/simong/EIC/detector_benchmarks_anl/sim_output/beamline/acceptanceTestXS2.edm4hep.root",
                         TString outFile            = "output.root",
                         std::string compactName    = "/home/simong/EIC/epic/install/share/epic/epic_ip6_extended.xml",
                         TString beampipeCanvasName = "acceptance_in_beampipe.png",

benchmarks/beamline/beamlineAnalysis.C

@@ -19,7 +19,7 @@
 using RVecS       = ROOT::VecOps::RVec<string>;
 using RNode       = ROOT::RDF::RNode;
 
-int beamlineAnalysis(   TString inFile          = "/scratch/EIC/G4out/beamline/beamlineTest.edm4hep.root",
+int beamlineAnalysis(   TString inFile          = "/home/simong/EIC/detector_benchmarks_anl/sim_output/beamline/acceptanceTestXS3.edm4hep.root",
                         TString outFile         = "output.root",
                         std::string compactName = "/home/simong/EIC/epic/install/share/epic/epic_ip6_extended.xml",
                         TString beamspotCanvasName = "beamspot_canvas.png",
@@ -79,6 +79,13 @@ int beamlineAnalysis(   TString inFile          = "/scratch/EIC/G4out/beamline/b
                 }
                 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"})
                 .Define("xdet",[](const ROOT::VecOps::RVec<volParams>& params) {
                 ROOT::VecOps::RVec<double> xPos;
                 for (const auto& param : params) {
@@ -170,6 +177,7 @@ int beamlineAnalysis(   TString inFile          = "/scratch/EIC/G4out/beamline/b
     std::map<TString,double> pipeXPos;
     std::map<TString,double> pipeZPos;
     std::map<TString,double> pipeRotation;
+    std::map<TString,bool> pipeIsConeSegment;
 
     // Queue up all actions
     auto xmin_ptr  = d1.Min("xpos");
@@ -201,6 +209,7 @@ int beamlineAnalysis(   TString inFile          = "/scratch/EIC/G4out/beamline/b
                           .Define("xmomf","xmom[pipeID=="+std::to_string(i)+"]")
                           .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("xdetf","xdet[pipeID=="+std::to_string(i)+"]")
                           .Define("zdetf","zdet[pipeID=="+std::to_string(i)+"]")
                           .Define("rotationf","rotation[pipeID=="+std::to_string(i)+"]");
@@ -273,6 +282,7 @@ int beamlineAnalysis(   TString inFile          = "/scratch/EIC/G4out/beamline/b
         pipeXPos[name]     = filterDF.Max("xdetf").GetValue();
         pipeZPos[name]     = filterDF.Max("zdetf").GetValue();
         pipeRotation[name] = filterDF.Max("rotationf").GetValue();
+        pipeIsConeSegment[name] = filterDF.Max("isConeSegmentf").GetValue();
 
         //Fit gaussian to the x, y, px and py distributions
         auto xhist = hHistsxy[name]->ProjectionX();
@@ -325,12 +335,15 @@ int beamlineAnalysis(   TString inFile          = "/scratch/EIC/G4out/beamline/b
         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 (pipeIsConeSegment[name] && pipeRadius > 0) {
+            TEllipse *circle = new TEllipse(0,0,pipeRadius);
+            circle->SetLineColor(kRed);
+            circle->SetLineWidth(2);
+            circle->SetFillStyle(0);
+            circle->Draw("same");
+        }
 
         // Add zoomed version in the top-right corner
         TPad *pad = new TPad("zoomPad", "Zoomed View", 0.65, 0.65, 1.0, 1.0);

benchmarks/beamline/shared_functions.h

@@ -4,9 +4,11 @@
 #include "edm4hep/SimTrackerHitCollection.h"
 #include "edm4hep/SimCalorimeterHitCollection.h"
 #include "DD4hep/VolumeManager.h"
+#include "DD4hep/DetElement.h"
 #include "TFile.h"
 
 using RVecHits = ROOT::VecOps::RVec<edm4hep::SimTrackerHitData>;
+using namespace dd4hep;
 
 //-----------------------------------------------------------------------------------------
 // Grab Component functor
@@ -108,6 +110,7 @@ struct volParams{
   double yPos;
   double zPos;
   double rotation;
+  bool isConeSegment;
 };
 
 // Functor to get the volume element parameters from the CellID
@@ -127,13 +130,20 @@ struct getVolumeParametersFromCellID {
       const Double_t*   rotationMatrix = transform.GetRotationMatrix(); // Compute rotation angle around the Y-axis
       double rotationAngleY = std::atan2(rotationMatrix[2], rotationMatrix[8]); // atan2(r13, r33)
       auto volume = detelement.volume();
-      dd4hep::ConeSegment cone = volume.solid();
+      dd4hep::Solid solid = volume.solid();
+      bool isCone = solid.isValid() && std::string(solid->IsA()->GetName()) == "TGeoConeSeg";
+      double radius = 0.0;
+      if (isCone) {
+        dd4hep::ConeSegment cone = solid;
+        radius = cone.rMax1();
+      }
       volParams params{
-        cone.rMax1(),
+        radius,
         translation[0],
         translation[1],
         translation[2],
-        rotationAngleY
+        rotationAngleY,
+        isCone
       };
       result.push_back(params);
     }
@@ -175,4 +185,13 @@ TH1F* CreateFittedHistogram(const std::string& histName,
   hist->SetMarkerColor(kRed);
 
   return hist;
+}
+
+void printHierarchy(const DetElement& de, int level = 0) {
+  std::string indent(level * 2, ' ');
+  std::cout << indent << "- " << de.name() << " (ID: " << de.id() << ")\n";
+
+  for (const auto& [childName, child] : de.children()) {
+    printHierarchy(child, level + 1);
+  }
 }