diff --git a/benchmarks/beamline/acceptanceAnalysis.C b/benchmarks/beamline/acceptanceAnalysis.C index 0a4e5c01..800ba024 100644 --- a/benchmarks/beamline/acceptanceAnalysis.C +++ b/benchmarks/beamline/acceptanceAnalysis.C @@ -19,7 +19,7 @@ using RVecS = ROOT::VecOps::RVec; 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", diff --git a/benchmarks/beamline/beamlineAnalysis.C b/benchmarks/beamline/beamlineAnalysis.C index 08fd7618..08ae6e33 100644 --- a/benchmarks/beamline/beamlineAnalysis.C +++ b/benchmarks/beamline/beamlineAnalysis.C @@ -19,7 +19,7 @@ using RVecS = ROOT::VecOps::RVec; 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& params) { + ROOT::VecOps::RVec cones; + for (const auto& param : params) { + cones.push_back(param.isConeSegment); + } + return cones; + }, {"pipeParameters"}) .Define("xdet",[](const ROOT::VecOps::RVec& params) { ROOT::VecOps::RVec xPos; for (const auto& param : params) { @@ -170,6 +177,7 @@ int beamlineAnalysis( TString inFile = "/scratch/EIC/G4out/beamline/b std::map pipeXPos; std::map pipeZPos; std::map pipeRotation; + std::map 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); diff --git a/benchmarks/beamline/shared_functions.h b/benchmarks/beamline/shared_functions.h index fede4b84..b0bf1e8e 100644 --- a/benchmarks/beamline/shared_functions.h +++ b/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; +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); + } } \ No newline at end of file