From 08586ff4ef31c02ce2ebcc6b6428df4eb794952f Mon Sep 17 00:00:00 2001 From: simonge Date: Mon, 16 Sep 2024 15:40:20 +0100 Subject: [PATCH 01/36] Added files --- Snakefile | 1 + benchmarks/beamline/Snakefile | 25 +++++++++++++++++ benchmarks/beamline/analysis.C | 50 ++++++++++++++++++++++++++++++++++ benchmarks/beamline/functors.h | 36 ++++++++++++++++++++++++ 4 files changed, 112 insertions(+) create mode 100644 benchmarks/beamline/Snakefile create mode 100644 benchmarks/beamline/analysis.C create mode 100644 benchmarks/beamline/functors.h diff --git a/Snakefile b/Snakefile index d645066a..c001f8b0 100644 --- a/Snakefile +++ b/Snakefile @@ -51,6 +51,7 @@ include: "benchmarks/insert_tau/Snakefile" include: "benchmarks/femc_electron/Snakefile" include: "benchmarks/femc_photon/Snakefile" include: "benchmarks/femc_pi0/Snakefile" +include: "benchmarks/beamline/Snakefile" use_s3 = config["remote_provider"].lower() == "s3" use_xrootd = config["remote_provider"].lower() == "xrootd" diff --git a/benchmarks/beamline/Snakefile b/benchmarks/beamline/Snakefile new file mode 100644 index 00000000..70eaf753 --- /dev/null +++ b/benchmarks/beamline/Snakefile @@ -0,0 +1,25 @@ +rule beamline_sim: + input: + inputfile="/scratch/EIC/Events/el_beam_18.hepmc", + output: + "/scratch/EIC/G4out/beamline/beamlineTest.edm4hep.root", + log: + "/scratch/EIC/G4out/beamline/beamlineTest.edm4hep.root.log", + params: + N_EVENTS=10000 + shell: + """ + exec npsim \ + --runType batch \ + --inputFiles {input.inputfile} \ + --random.seed 1234 \ + --numberOfEvents {params.N_EVENTS} \ + --compactFile $DETECTOR_PATH/epic_ip6_extended.xml \ + --outputFile {output} + """ +rule beamline_analysis + input: + "/scratch/EIC/G4out/beamline/beamlineTest.edm4hep.root", + output: + "/scratch/EIC/ReconOut/beamline/beamlineTestAnalysis.root" + log: \ No newline at end of file diff --git a/benchmarks/beamline/analysis.C b/benchmarks/beamline/analysis.C new file mode 100644 index 00000000..a0551139 --- /dev/null +++ b/benchmarks/beamline/analysis.C @@ -0,0 +1,50 @@ +// Script to plot the x and y positions and momentum of beamline particles as they pass through the magnets + +#include +#include "DD4hep/Detector.h" +#include "DDRec/CellIDPositionConverter.h" +#include "edm4hep/MCParticleCollection.h" +#include "edm4hep/SimTrackerHit.h" +#include "edm4hep/SimCalorimeterHit.h" +#include "ROOT/RDataFrame.hxx" +#include "ROOT/RDF/RInterface.hxx" +#include "ROOT/RVec.hxx" +#include "functors.h" + +using RVecS = ROOT::VecOps::RVec; +using RNode = ROOT::RDF::RNode; + +void analysis( std::string inFile = "/scratch/EIC/G4out/beamline/beamlineTest.edm4hep.root", + std::string compactName = "/home/simong/EIC/epic/install/share/epic/epic_ip6_extended.xml"){ + + //Set implicit multi-threading + ROOT::EnableImplicitMT(); + + //Load the detector config + dd4hep::Detector& detector = dd4hep::Detector::getInstance(); + detector.fromCompact(compactName); + + ROOT::RDataFrame d0("events",inFile, {"BackwardsBeamlineHits"}); + RNode d1 = d0; + RVecS colNames = d1.GetColumnNames(); + + //Get the collection + std::string readoutName = "BackwardsBeamlineHits"; + + if(Any(colNames==readoutName)){ + + auto ids = detector.readout(readoutName).idSpec().fields(); + + for(auto &[key,id]: ids){ + TString colName = key+"ID"; + d1 = d1.Define(colName,getSubID(key,detector),{readoutName}); + } + } + else{ + std::cout << "Collection " << readoutName << " not found in file" << std::endl; + return; + } + + d1.Snapshot("events","output.root"); + +} diff --git a/benchmarks/beamline/functors.h b/benchmarks/beamline/functors.h new file mode 100644 index 00000000..d81b4567 --- /dev/null +++ b/benchmarks/beamline/functors.h @@ -0,0 +1,36 @@ +#include "DD4hep/Detector.h" +#include "DDRec/CellIDPositionConverter.h" +#include "edm4hep/MCParticleCollection.h" +#include "edm4hep/SimTrackerHit.h" +#include "edm4hep/SimCalorimeterHit.h" + +using RVecHits = ROOT::VecOps::RVec; + +//----------------------------------------------------------------------------------------- +// Grab Component functor +//----------------------------------------------------------------------------------------- +struct getSubID{ + getSubID(std::string cname, dd4hep::Detector& det, std::string rname = "BackwardsBeamlineHits") : componentName(cname), detector(det), readoutName(rname){} + + ROOT::VecOps::RVec operator()(const RVecHits& evt) { + auto decoder = detector.readout(readoutName).idSpec().decoder(); + auto indexID = decoder->index(componentName); + ROOT::VecOps::RVec result; + for(const auto& h: evt) { + result.push_back(decoder->get(h.cellID,indexID)); + } + return result; + }; + + void SetComponent(std::string cname){ + componentName = cname; + } + void SetReadout(std::string rname){ + readoutName = rname; + } + + private: + std::string componentName; + dd4hep::Detector& detector; + std::string readoutName; +}; \ No newline at end of file From df1b38ee23b0c43909c72e3fc5976ce40f68740a Mon Sep 17 00:00:00 2001 From: simonge Date: Wed, 30 Oct 2024 10:18:21 +0000 Subject: [PATCH 02/36] Update --- benchmarks/beamline/Snakefile | 25 ++++--- benchmarks/beamline/analysis.C | 116 +++++++++++++++++++++++++++++---- benchmarks/beamline/functors.h | 71 +++++++++++++++++++- 3 files changed, 189 insertions(+), 23 deletions(-) diff --git a/benchmarks/beamline/Snakefile b/benchmarks/beamline/Snakefile index 70eaf753..94b5c1e7 100644 --- a/benchmarks/beamline/Snakefile +++ b/benchmarks/beamline/Snakefile @@ -1,25 +1,30 @@ rule beamline_sim: input: - inputfile="/scratch/EIC/Events/el_beam_18.hepmc", + macro="beamGPS.mac", output: "/scratch/EIC/G4out/beamline/beamlineTest.edm4hep.root", log: "/scratch/EIC/G4out/beamline/beamlineTest.edm4hep.root.log", - params: - N_EVENTS=10000 shell: """ exec npsim \ - --runType batch \ - --inputFiles {input.inputfile} \ - --random.seed 1234 \ - --numberOfEvents {params.N_EVENTS} \ + --runType run \ + --enableG4GPS \ + --macroFile {input.macro} \ --compactFile $DETECTOR_PATH/epic_ip6_extended.xml \ - --outputFile {output} + --outputFile {output} \ """ -rule beamline_analysis + +rule beamline_analysis: input: "/scratch/EIC/G4out/beamline/beamlineTest.edm4hep.root", output: "/scratch/EIC/ReconOut/beamline/beamlineTestAnalysis.root" - log: \ No newline at end of file + log: + "/scratch/EIC/ReconOut/beamline/beamlineTestAnalysis.root.log" + params: + xml=os.getenv("DETECTOR_PATH")+"/epic_ip6_extended.xml", + shell: + """ + root -l -b -q 'analysis.C("{input}", "{output}", "{params.xml}")' + """ diff --git a/benchmarks/beamline/analysis.C b/benchmarks/beamline/analysis.C index a0551139..0c4c9465 100644 --- a/benchmarks/beamline/analysis.C +++ b/benchmarks/beamline/analysis.C @@ -4,22 +4,35 @@ #include "DD4hep/Detector.h" #include "DDRec/CellIDPositionConverter.h" #include "edm4hep/MCParticleCollection.h" -#include "edm4hep/SimTrackerHit.h" -#include "edm4hep/SimCalorimeterHit.h" +#include "edm4hep/SimTrackerHitCollection.h" +#include "edm4hep/SimCalorimeterHitCollection.h" +#include "edm4hep/Vector3f.h" +#include "edm4hep/Vector3d.h" #include "ROOT/RDataFrame.hxx" #include "ROOT/RDF/RInterface.hxx" #include "ROOT/RVec.hxx" #include "functors.h" +#include "TCanvas.h" +#include "TStyle.h" using RVecS = ROOT::VecOps::RVec; using RNode = ROOT::RDF::RNode; -void analysis( std::string inFile = "/scratch/EIC/G4out/beamline/beamlineTest.edm4hep.root", +void analysis( TString inFile = "/scratch/EIC/G4out/beamline/beamlineTest.edm4hep.root", + TString outFile = "output.root", std::string compactName = "/home/simong/EIC/epic/install/share/epic/epic_ip6_extended.xml"){ + //Set ROOT style + gStyle->SetPadLeftMargin(0.25); // Set left margin + gStyle->SetPadRightMargin(0.15); // Set right margin + gStyle->SetPadTopMargin(0.05); // Set top margin + gStyle->SetPadBottomMargin(0.15);// Set bottom margin + gStyle->SetTitleYOffset(2); // Adjust y-axis title offset + gStyle->SetOptStat(1110); + //Set implicit multi-threading ROOT::EnableImplicitMT(); - + //Load the detector config dd4hep::Detector& detector = dd4hep::Detector::getInstance(); detector.fromCompact(compactName); @@ -27,24 +40,105 @@ void analysis( std::string inFile = "/scratch/EIC/G4out/beamline/beamlineT ROOT::RDataFrame d0("events",inFile, {"BackwardsBeamlineHits"}); RNode d1 = d0; RVecS colNames = d1.GetColumnNames(); + + //Set number of entries to process + // d1 = d1.Range(0,1000); //Get the collection std::string readoutName = "BackwardsBeamlineHits"; + std::cout << "Running lazy RDataframe execution" << std::endl; + if(Any(colNames==readoutName)){ - auto ids = detector.readout(readoutName).idSpec().fields(); + d1 = d1.Define("pipeID",getSubID("pipe",detector),{readoutName}) + .Define("endID",getSubID("end",detector),{readoutName}); + + //global x,y,z position and momentum + d1 = d1.Define("xpos_global","BackwardsBeamlineHits.position.x") + .Define("ypos_global","BackwardsBeamlineHits.position.y") + .Define("zpos_global","BackwardsBeamlineHits.position.z") + .Define("px_global","BackwardsBeamlineHits.momentum.x") + .Define("py_global","BackwardsBeamlineHits.momentum.y") + .Define("pz_global","BackwardsBeamlineHits.momentum.z"); + + d1 = d1.Define("hitPosMom",globalToLocal(detector),{readoutName}) + .Define("xpos","hitPosMom[0]") + .Define("ypos","hitPosMom[1]") + .Define("zpos","hitPosMom[2]") + .Define("xmom","hitPosMom[3]") + .Define("ymom","hitPosMom[4]") + .Define("zmom","hitPosMom[5]"); - for(auto &[key,id]: ids){ - TString colName = key+"ID"; - d1 = d1.Define(colName,getSubID(key,detector),{readoutName}); - } } else{ std::cout << "Collection " << readoutName << " not found in file" << std::endl; return; + } + + std::cout << "Executing Analysis and creating histograms" << std::endl; + + //Create array of histogram results + std::map> hHistsxpx; + std::map> hHistsypy; + + + //Create histograms + for(int i=0; i<=7; i++){ + + std::string name = "pipeID"; + name += std::to_string(i); + auto filterDF = d1.Define("xposf","xpos["+std::to_string(i)+"]") + .Define("yposf","ypos["+std::to_string(i)+"]") + .Define("xmomf","xmom["+std::to_string(i)+"]") + .Define("ymomf","ymom["+std::to_string(i)+"]"); + + //Calculate Min and Max values + auto xmin = filterDF.Min("xposf").GetValue(); + auto xmax = filterDF.Max("xposf").GetValue(); + auto ymin = filterDF.Min("yposf").GetValue(); + auto ymax = filterDF.Max("yposf").GetValue(); + auto pxmin = filterDF.Min("xmomf").GetValue(); + auto pxmax = filterDF.Max("xmomf").GetValue(); + auto pymin = filterDF.Min("ymomf").GetValue(); + auto pymax = filterDF.Max("ymomf").GetValue(); + + TString xname = name+";x offset [mm]; x trajectory component"; + TString yname = name+";y offset [mm]; y trajectory component"; + hHistsxpx[name] = filterDF.Histo2D({xname,xname,100,xmin,xmax,100,pxmin,pxmax},"xposf","xmomf"); + hHistsypy[name] = filterDF.Histo2D({yname,yname,100,ymin,ymax,100,pymin,pymax},"yposf","ymomf"); + } - - d1.Snapshot("events","output.root"); + TCanvas *cX = new TCanvas("c1","c1",3000,1600); + cX->Divide(4,2); + + int i=1; + //Write histograms to file + for(auto& h : hHistsxpx){ + cX->cd(i++); + h.second->Draw("colz"); + } + + + TCanvas *cY = new TCanvas("c2","c2",3000,1600); + cY->Divide(4,2); + + i=1; + for(auto& h : hHistsypy){ + cY->cd(i++); + h.second->Draw("colz"); + } + + TFile *f = new TFile(outFile,"RECREATE"); + cX->Write(); + cY->Write(); + + f->Close(); + + std::cout << "Saving events to file" << std::endl; + + // ROOT::RDF::RSnapshotOptions opts; + // opts.fMode = "UPDATE"; + // d1.Snapshot("events",outFile,{"pipeID","endID","xpos","ypos","zpos","xmom","ymom","zmom","xpos_global","ypos_global","zpos_global","px_global","py_global","pz_global"},opts); } diff --git a/benchmarks/beamline/functors.h b/benchmarks/beamline/functors.h index d81b4567..ee2414d0 100644 --- a/benchmarks/beamline/functors.h +++ b/benchmarks/beamline/functors.h @@ -1,8 +1,10 @@ #include "DD4hep/Detector.h" #include "DDRec/CellIDPositionConverter.h" #include "edm4hep/MCParticleCollection.h" -#include "edm4hep/SimTrackerHit.h" -#include "edm4hep/SimCalorimeterHit.h" +#include "edm4hep/SimTrackerHitCollection.h" +#include "edm4hep/SimCalorimeterHitCollection.h" +#include "DD4hep/VolumeManager.h" +#include "TFile.h" using RVecHits = ROOT::VecOps::RVec; @@ -33,4 +35,69 @@ struct getSubID{ std::string componentName; dd4hep::Detector& detector; std::string readoutName; +}; + +//----------------------------------------------------------------------------------------- +// Transform global x,y,z position and momentum into local coordinates +//----------------------------------------------------------------------------------------- +struct globalToLocal{ + globalToLocal(dd4hep::Detector& det) : detector(det){ + volumeManager = dd4hep::VolumeManager::getVolumeManager(det); + } + + ROOT::VecOps::RVec> operator()(const RVecHits& evt) { + + ROOT::VecOps::RVec> result; + ROOT::VecOps::RVec xvec; + ROOT::VecOps::RVec yvec; + ROOT::VecOps::RVec zvec; + ROOT::VecOps::RVec pxvec; + ROOT::VecOps::RVec pyvec; + ROOT::VecOps::RVec pzvec; + + for(const auto& h: evt) { + auto cellID = h.cellID; + // dd4hep::DetElement detelement = volumeManager.lookupDetElement(cellID); + // auto detelement = volumeManager.lookupVolumePlacement(cellID); + auto detelement = volumeManager.lookupDetElement(cellID); + const TGeoMatrix& transform = detelement.nominal().worldTransformation(); + // transform.Print(); + //auto transform = volumeManager.worldTransformation(cellID); + //Convert position to local coordinates + auto pos = h.position; + Double_t globalPos[3] = {pos.x/10, pos.y/10, pos.z/10}; + Double_t localPos[3]; + transform.MasterToLocal(globalPos, localPos); + // std::cout << "Global: " << globalPos[0] << " " << globalPos[1] << " " << globalPos[2] << std::endl; + // std::cout << "Local: " << localPos[0] << " " << localPos[1] << " " << localPos[2] << std::endl; + //Transform global momentum to local coordinates + auto mom = h.momentum; + Double_t globalMom[3] = {mom.x, mom.y, mom.z}; + Double_t localMom[3]; + transform.MasterToLocalVect(globalMom, localMom); + // std::cout << "Global: " << globalMom[0] << " " << globalMom[1] << " " << globalMom[2] << std::endl; + // std::cout << "Local: " << localMom[0] << " " << localMom[1] << " " << localMom[2] << std::endl; + + xvec.push_back(localPos[0]); + yvec.push_back(localPos[1]); + zvec.push_back(localPos[2]); + pxvec.push_back(localMom[0]); + pyvec.push_back(localMom[1]); + pzvec.push_back(localMom[2]); + + } + + result.push_back(xvec); + result.push_back(yvec); + result.push_back(zvec); + result.push_back(pxvec); + result.push_back(pyvec); + result.push_back(pzvec); + + return result; + }; + + private: + dd4hep::Detector& detector; + dd4hep::VolumeManager volumeManager; }; \ No newline at end of file From 810cc3b84ee3ff440de661fe9fd87e5672766112 Mon Sep 17 00:00:00 2001 From: simonge Date: Tue, 11 Mar 2025 17:42:35 +0000 Subject: [PATCH 03/36] Make canvas name more useful --- benchmarks/beamline/Snakefile | 10 +++++----- benchmarks/beamline/analysis.C | 4 ++-- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/benchmarks/beamline/Snakefile b/benchmarks/beamline/Snakefile index 94b5c1e7..13cd66a7 100644 --- a/benchmarks/beamline/Snakefile +++ b/benchmarks/beamline/Snakefile @@ -2,9 +2,9 @@ rule beamline_sim: input: macro="beamGPS.mac", output: - "/scratch/EIC/G4out/beamline/beamlineTest.edm4hep.root", + "/scratch/EIC/G4out/beamline/beamlineTest{tag}.edm4hep.root", log: - "/scratch/EIC/G4out/beamline/beamlineTest.edm4hep.root.log", + "/scratch/EIC/G4out/beamline/beamlineTest{tag}.edm4hep.root.log", shell: """ exec npsim \ @@ -17,11 +17,11 @@ rule beamline_sim: rule beamline_analysis: input: - "/scratch/EIC/G4out/beamline/beamlineTest.edm4hep.root", + "/scratch/EIC/G4out/beamline/beamlineTest{tag}.edm4hep.root", output: - "/scratch/EIC/ReconOut/beamline/beamlineTestAnalysis.root" + "/scratch/EIC/ReconOut/beamline/beamlineTestAnalysis{tag}.root" log: - "/scratch/EIC/ReconOut/beamline/beamlineTestAnalysis.root.log" + "/scratch/EIC/ReconOut/beamline/beamlineTestAnalysis{tag}.root.log" params: xml=os.getenv("DETECTOR_PATH")+"/epic_ip6_extended.xml", shell: diff --git a/benchmarks/beamline/analysis.C b/benchmarks/beamline/analysis.C index 0c4c9465..77672054 100644 --- a/benchmarks/beamline/analysis.C +++ b/benchmarks/beamline/analysis.C @@ -110,7 +110,7 @@ void analysis( TString inFile = "/scratch/EIC/G4out/beamline/beamlineTest. } - TCanvas *cX = new TCanvas("c1","c1",3000,1600); + TCanvas *cX = new TCanvas("x_px_canvas","x_px_canvas",3000,1600); cX->Divide(4,2); int i=1; @@ -121,7 +121,7 @@ void analysis( TString inFile = "/scratch/EIC/G4out/beamline/beamlineTest. } - TCanvas *cY = new TCanvas("c2","c2",3000,1600); + TCanvas *cY = new TCanvas("y_py_canvas","y_py_canvas",3000,1600); cY->Divide(4,2); i=1; From 0973e5d250c43e82b6bf2ddb508688485af3e1c5 Mon Sep 17 00:00:00 2001 From: simonge Date: Tue, 11 Mar 2025 17:54:21 +0000 Subject: [PATCH 04/36] Add beam generator macro --- benchmarks/beamline/beamGPS.mac | 10 ++++++++++ 1 file changed, 10 insertions(+) create mode 100644 benchmarks/beamline/beamGPS.mac diff --git a/benchmarks/beamline/beamGPS.mac b/benchmarks/beamline/beamGPS.mac new file mode 100644 index 00000000..28d771ba --- /dev/null +++ b/benchmarks/beamline/beamGPS.mac @@ -0,0 +1,10 @@ +/gps/ang/type beam2d +/gps/ang/sigma_x 0.0002017 rad # 201.7 urad +/gps/ang/sigma_y 0.0001873 rad # 187.3 urad +/gps/pos/type Beam +/gps/pos/sigma_x 0.119 mm +/gps/pos/sigma_y 0.0107 mm +/gps/energy 18 GeV +/gps/particle e- + +/run/beamOn 100000 \ No newline at end of file From 1323fe045148f9142e3e90d448794b6da7e9a5f3 Mon Sep 17 00:00:00 2001 From: simonge Date: Wed, 30 Apr 2025 18:26:15 +0100 Subject: [PATCH 05/36] Update analysis --- benchmarks/beamline/analysis.C | 373 +++++++++++++++++++++++++++++--- benchmarks/beamline/beamGPS.mac | 3 +- benchmarks/beamline/functors.h | 77 ++++++- 3 files changed, 419 insertions(+), 34 deletions(-) diff --git a/benchmarks/beamline/analysis.C b/benchmarks/beamline/analysis.C index 77672054..f37c2ce1 100644 --- a/benchmarks/beamline/analysis.C +++ b/benchmarks/beamline/analysis.C @@ -15,6 +15,7 @@ #include "TCanvas.h" #include "TStyle.h" + using RVecS = ROOT::VecOps::RVec; using RNode = ROOT::RDF::RNode; @@ -23,12 +24,17 @@ void analysis( TString inFile = "/scratch/EIC/G4out/beamline/beamlineTest. std::string compactName = "/home/simong/EIC/epic/install/share/epic/epic_ip6_extended.xml"){ //Set ROOT style - gStyle->SetPadLeftMargin(0.25); // Set left margin - gStyle->SetPadRightMargin(0.15); // Set right margin - gStyle->SetPadTopMargin(0.05); // Set top margin - gStyle->SetPadBottomMargin(0.15);// Set bottom margin - gStyle->SetTitleYOffset(2); // Adjust y-axis title offset - gStyle->SetOptStat(1110); + gStyle->SetPadLeftMargin(0.1); // Set left margin + gStyle->SetPadRightMargin(0.0); // Set right margin + gStyle->SetPadTopMargin(0.0); // Set top margin + gStyle->SetPadBottomMargin(0.1);// Set bottom margin + gStyle->SetTitleAlign(13); + gStyle->SetTitleX(0.12); // Place the title on the top right + gStyle->SetTitleY(0.985); // Place the title on the top right + gStyle->SetTitleSize(0.08, "t"); + gStyle->SetTitleXOffset(1.0); // Adjust y-axis title offset + gStyle->SetTitleYOffset(1.0); // Adjust y-axis title offset + gStyle->SetOptStat(0); //Set implicit multi-threading ROOT::EnableImplicitMT(); @@ -45,14 +51,44 @@ void analysis( TString inFile = "/scratch/EIC/G4out/beamline/beamlineTest. // d1 = d1.Range(0,1000); //Get the collection - std::string readoutName = "BackwardsBeamlineHits"; + std::string readoutName = "BackwardsBeamlineHits"; std::cout << "Running lazy RDataframe execution" << std::endl; if(Any(colNames==readoutName)){ d1 = d1.Define("pipeID",getSubID("pipe",detector),{readoutName}) - .Define("endID",getSubID("end",detector),{readoutName}); + .Define("endID",getSubID("end",detector),{readoutName}) + .Define("pipeParameters",getVolumeParametersFromCellID(detector),{readoutName}) + .Define("pipeRadius",[](const ROOT::VecOps::RVec& params) { + ROOT::VecOps::RVec radii; + for (const auto& param : params) { + radii.push_back(param.radius); + } + return radii; + }, {"pipeParameters"}) + .Define("xdet",[](const ROOT::VecOps::RVec& params) { + ROOT::VecOps::RVec xPos; + for (const auto& param : params) { + xPos.push_back(param.xPos); + } + return xPos; + }, {"pipeParameters"}) + .Define("zdet",[](const ROOT::VecOps::RVec& params) { + ROOT::VecOps::RVec zPos; + for (const auto& param : params) { + zPos.push_back(param.zPos); + } + return zPos; + }, {"pipeParameters"}) + .Define("rotation",[](const ROOT::VecOps::RVec& params) { + ROOT::VecOps::RVec rotation; + for (const auto& param : params) { + rotation.push_back(param.rotation); + } + return rotation; + }, {"pipeParameters"}); + //global x,y,z position and momentum d1 = d1.Define("xpos_global","BackwardsBeamlineHits.position.x") @@ -66,9 +102,13 @@ void analysis( TString inFile = "/scratch/EIC/G4out/beamline/beamlineTest. .Define("xpos","hitPosMom[0]") .Define("ypos","hitPosMom[1]") .Define("zpos","hitPosMom[2]") - .Define("xmom","hitPosMom[3]") - .Define("ymom","hitPosMom[4]") - .Define("zmom","hitPosMom[5]"); + .Define("xmomMag","hitPosMom[3]") + .Define("ymomMag","hitPosMom[4]") + .Define("zmomMag","hitPosMom[5]") + .Define("momMag","sqrt(xmomMag*xmomMag+ymomMag*ymomMag+zmomMag*zmomMag)") + .Define("xmom","xmomMag/momMag") + .Define("ymom","ymomMag/momMag") + .Define("zmom","zmomMag/momMag"); } else{ @@ -76,13 +116,58 @@ void analysis( TString inFile = "/scratch/EIC/G4out/beamline/beamlineTest. return; } + // Calculate the maximum pipe radius for plotting + auto maxPipeRadius = 1.2*d1.Max("pipeRadius").GetValue(); + std::cout << "Executing Analysis and creating histograms" << std::endl; //Create array of histogram results + std::map> hHistsxy; + std::map> hHistsxyZoom; std::map> hHistsxpx; std::map> hHistsypy; + std::map xMeans; + std::map yMeans; + std::map xStdDevs; + std::map yStdDevs; + std::map pxMeans; + std::map pyMeans; + std::map pxStdDevs; + std::map pyStdDevs; + + //Fit paremeter and error maps + std::map xMeanFit; + std::map yMeanFit; + std::map xMeanFitErr; + std::map yMeanFitErr; + std::map xStdDevFit; + std::map yStdDevFit; + std::map xStdDevFitErr; + std::map yStdDevFitErr; + std::map pxMeanFit; + std::map pyMeanFit; + std::map pxMeanFitErr; + std::map pyMeanFitErr; + std::map pxStdDevFit; + std::map pyStdDevFit; + std::map pxStdDevFitErr; + std::map pyStdDevFitErr; + + std::map pipeRadii; + std::map pipeXPos; + std::map pipeZPos; + std::map pipeRotation; + auto xmin = d1.Min("xpos").GetValue(); + auto xmax = d1.Max("xpos").GetValue(); + auto ymin = d1.Min("ypos").GetValue(); + auto ymax = d1.Max("ypos").GetValue(); + auto pxmin = d1.Min("xmom").GetValue(); + auto pxmax = d1.Max("xmom").GetValue(); + auto pymin = d1.Min("ymom").GetValue(); + auto pymax = d1.Max("ymom").GetValue(); + //Create histograms for(int i=0; i<=7; i++){ @@ -91,54 +176,278 @@ void analysis( TString inFile = "/scratch/EIC/G4out/beamline/beamlineTest. auto filterDF = d1.Define("xposf","xpos["+std::to_string(i)+"]") .Define("yposf","ypos["+std::to_string(i)+"]") .Define("xmomf","xmom["+std::to_string(i)+"]") - .Define("ymomf","ymom["+std::to_string(i)+"]"); + .Define("ymomf","ymom["+std::to_string(i)+"]") + .Define("pipeRadiusf","pipeRadius["+std::to_string(i)+"]") + .Define("xdetf","xdet["+std::to_string(i)+"]") + .Define("zdetf","zdet["+std::to_string(i)+"]") + .Define("rotationf","rotation["+std::to_string(i)+"]"); + //Calculate Min and Max values - auto xmin = filterDF.Min("xposf").GetValue(); - auto xmax = filterDF.Max("xposf").GetValue(); - auto ymin = filterDF.Min("yposf").GetValue(); - auto ymax = filterDF.Max("yposf").GetValue(); - auto pxmin = filterDF.Min("xmomf").GetValue(); - auto pxmax = filterDF.Max("xmomf").GetValue(); - auto pymin = filterDF.Min("ymomf").GetValue(); - auto pymax = filterDF.Max("ymomf").GetValue(); - - TString xname = name+";x offset [mm]; x trajectory component"; - TString yname = name+";y offset [mm]; y trajectory component"; - hHistsxpx[name] = filterDF.Histo2D({xname,xname,100,xmin,xmax,100,pxmin,pxmax},"xposf","xmomf"); - hHistsypy[name] = filterDF.Histo2D({yname,yname,100,ymin,ymax,100,pymin,pymax},"yposf","ymomf"); + auto xminf = filterDF.Min("xposf").GetValue(); + auto xmaxf = filterDF.Max("xposf").GetValue(); + auto yminf = filterDF.Min("yposf").GetValue(); + auto ymaxf = filterDF.Max("yposf").GetValue(); + auto pxminf = filterDF.Min("xmomf").GetValue(); + auto pxmaxf = filterDF.Max("xmomf").GetValue(); + auto pyminf = filterDF.Min("ymomf").GetValue(); + auto pymaxf = filterDF.Max("ymomf").GetValue(); + // Calculate means and standard deviations + xMeans[name] = filterDF.Mean("xposf").GetValue(); + yMeans[name] = filterDF.Mean("yposf").GetValue(); + xStdDevs[name] = filterDF.StdDev("xposf").GetValue(); + yStdDevs[name] = filterDF.StdDev("yposf").GetValue(); + pxMeans[name] = filterDF.Mean("xmomf").GetValue(); + pyMeans[name] = filterDF.Mean("ymomf").GetValue(); + pxStdDevs[name] = filterDF.StdDev("xmomf").GetValue(); + pyStdDevs[name] = filterDF.StdDev("ymomf").GetValue(); + + // Calculate axes for zoomed beamspot histogram so that it is quare around the mean x and y + double halfrange = std::max({xMeans[name]-xminf, xmaxf-xMeans[name], yMeans[name]-yminf, ymaxf-yMeans[name]}); + double xMinZoom = xMeans[name] - halfrange; + double xMaxZoom = xMeans[name] + halfrange; + double yMinZoom = yMeans[name] - halfrange; + double yMaxZoom = yMeans[name] + halfrange; + + TString beamspotName = "Beamspot ID"+std::to_string(i)+";x offset [cm]; y offset [cm]"; + TString xyname = name+";x Offset [cm]; y Offset [cm]"; + TString xname = name+";x Offset [cm]; x trajectory component"; + TString yname = name+";y Offset [cm]; y trajectory component"; + hHistsxy[name] = filterDF.Histo2D({beamspotName,beamspotName,400,-maxPipeRadius,maxPipeRadius,400,-maxPipeRadius,maxPipeRadius},"xposf","yposf"); + hHistsxyZoom[name] = filterDF.Histo2D({xyname,xyname,100,xMinZoom,xMaxZoom,100,yMinZoom,yMaxZoom},"xposf","yposf"); + hHistsxpx[name] = filterDF.Histo2D({xname,xname,400,xmin,xmax,400,pxmin,pxmax},"xposf","xmomf"); + hHistsypy[name] = filterDF.Histo2D({yname,yname,400,ymin,ymax,400,pymin,pymax},"yposf","ymomf"); + + //Parameters of the pipe + pipeRadii[name] = filterDF.Max("pipeRadiusf").GetValue(); + pipeXPos[name] = filterDF.Max("xdetf").GetValue(); + pipeZPos[name] = filterDF.Max("zdetf").GetValue(); + pipeRotation[name] = filterDF.Max("rotationf").GetValue(); + //Fit gaussian to the x, y, px and py distributions + auto xhist = hHistsxy[name]->ProjectionX(); + auto yhist = hHistsxy[name]->ProjectionY(); + auto pxhist = hHistsxpx[name]->ProjectionY(); + auto pyhist = hHistsypy[name]->ProjectionY(); + xhist->Fit("gaus","Q"); + yhist->Fit("gaus","Q"); + pxhist->Fit("gaus","Q"); + pyhist->Fit("gaus","Q"); + //Get the fit parameters and errors + xMeanFit[name] = xhist->GetFunction("gaus")->GetParameter(1); + yMeanFit[name] = yhist->GetFunction("gaus")->GetParameter(1); + xMeanFitErr[name] = xhist->GetFunction("gaus")->GetParError(1); + yMeanFitErr[name] = yhist->GetFunction("gaus")->GetParError(1); + xStdDevFit[name] = xhist->GetFunction("gaus")->GetParameter(2); + yStdDevFit[name] = yhist->GetFunction("gaus")->GetParameter(2); + xStdDevFitErr[name] = xhist->GetFunction("gaus")->GetParError(2); + yStdDevFitErr[name] = yhist->GetFunction("gaus")->GetParError(2); + pxMeanFit[name] = pxhist->GetFunction("gaus")->GetParameter(1); + pyMeanFit[name] = pyhist->GetFunction("gaus")->GetParameter(1); + pxMeanFitErr[name] = pxhist->GetFunction("gaus")->GetParError(1); + pyMeanFitErr[name] = pyhist->GetFunction("gaus")->GetParError(1); + pxStdDevFit[name] = pxhist->GetFunction("gaus")->GetParameter(2); + pyStdDevFit[name] = pyhist->GetFunction("gaus")->GetParameter(2); + pxStdDevFitErr[name] = pxhist->GetFunction("gaus")->GetParError(2); + pyStdDevFitErr[name] = pyhist->GetFunction("gaus")->GetParError(2); + } + // Create histograms of the beamspot + TCanvas *cXY = new TCanvas("beamspot_canvas","beamspot_canvas",3000,1600); + cXY->Divide(4,2); + int i=1; + for(auto [name,h] : hHistsxy){ + // Get the pipe radius for this histogram + auto pipeRadius = pipeRadii[name]; + 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"); + + // Add zoomed version in the top-right corner + TPad *pad = new TPad("zoomPad", "Zoomed View", 0.65, 0.65, 1.0, 1.0); + pad->SetFillStyle(0); // Transparent background + pad->Draw(); + pad->cd(); + // Draw the zoomed histogram without its title or axis titles + hHistsxyZoom[name]->SetTitle(""); + hHistsxyZoom[name]->GetXaxis()->SetTitle(""); + hHistsxyZoom[name]->GetYaxis()->SetTitle(""); + hHistsxyZoom[name]->Draw("col"); + cXY->cd(); // Return to the main canvas + } + + // x-px canvas TCanvas *cX = new TCanvas("x_px_canvas","x_px_canvas",3000,1600); cX->Divide(4,2); - int i=1; + i=1; //Write histograms to file for(auto& h : hHistsxpx){ cX->cd(i++); - h.second->Draw("colz"); + h.second->Draw("col"); } - + // y-py canvas TCanvas *cY = new TCanvas("y_py_canvas","y_py_canvas",3000,1600); cY->Divide(4,2); i=1; for(auto& h : hHistsypy){ cY->cd(i++); - h.second->Draw("colz"); + h.second->Draw("col"); } + // Save 2D canvases as pngs + cXY->SaveAs("beamspot.png"); + cX->SaveAs("x_px.png"); + cY->SaveAs("y_py.png"); + + // --------------------------------------------------------------------------- + // Create histograms showing the fitted means and standard deviations of the positions and momenta + // --------------------------------------------------------------------------- + + // Create histograms for fitted X means and standard deviations + TH1F* hFittedXMeans = CreateFittedHistogram("hFittedXMeans", + "Mean X Offset [cm]", + xMeanFit, + xMeanFitErr, + "Pipe ID"); + + TH1F* hFittedXStdDevs = CreateFittedHistogram("hFittedXStdDevs", + "Std Deviation X Offset [cm]", + xStdDevFit, + xStdDevFitErr, + "Pipe ID"); + + // Create histograms for fitted Y means and standard deviations + TH1F* hFittedYMeans = CreateFittedHistogram("hFittedYMeans", + "Mean Y Offset [cm]", + yMeanFit, + yMeanFitErr, + "Pipe ID"); + + TH1F* hFittedYStdDevs = CreateFittedHistogram("hFittedYStdDevs", + "Std Deviation Y Offset [cm]", + yStdDevFit, + yStdDevFitErr, + "Pipe ID"); + + TH1F* hFittedPxMeans = CreateFittedHistogram("hFittedPxMeans", + "Mean Px", + pxMeanFit, + pxMeanFitErr, + "Pipe ID"); + TH1F* hFittedPyMeans = CreateFittedHistogram("hFittedPyMeans", + "Mean Py", + pyMeanFit, + pyMeanFitErr, + "Pipe ID"); + TH1F* hFittedPxStdDevs = CreateFittedHistogram("hFittedPxStdDevs", + "Std Deviation Px", + pxStdDevFit, + pxStdDevFitErr, + "Pipe ID"); + TH1F* hFittedPyStdDevs = CreateFittedHistogram("hFittedPyStdDevs", + "Std Deviation Py", + pyStdDevFit, + pyStdDevFitErr, + "Pipe ID"); + + + // Create a canvas for the fitted beamspot means and standard deviations + TCanvas *cFittedMeans = new TCanvas("cFittedMeans", "Fitted Beamspot Means and Std Deviation", 1200, 800); + cFittedMeans->Divide(2, 2); + cFittedMeans->cd(1); + hFittedXMeans->Draw("E1"); // "E1" draws error bars + cFittedMeans->cd(2); + hFittedXStdDevs->Draw("E1"); // "E1" draws error bars + cFittedMeans->cd(3); + hFittedYMeans->Draw("E1"); // "E1" draws error bars + cFittedMeans->cd(4); + hFittedYStdDevs->Draw("E1"); // "E1" draws error bars + cFittedMeans->SetGrid(); + cFittedMeans->Update(); + // Save the canvas as a PNG file + cFittedMeans->SaveAs("fitted_means_stddevs.png"); + + // Create a canvas for the fitted momentum means and standard deviations + TCanvas *cFittedMomentumMeans = new TCanvas("cFittedMomentumMeans", "Fitted Momentum Means and Std Deviation", 1200, 800); + cFittedMomentumMeans->Divide(2, 2); + cFittedMomentumMeans->cd(1); + hFittedPxMeans->Draw("E1"); // "E1" draws error bars + cFittedMomentumMeans->cd(2); + hFittedPxStdDevs->Draw("E1"); // "E1" draws error bars + cFittedMomentumMeans->cd(3); + hFittedPyMeans->Draw("E1"); // "E1" draws error bars + cFittedMomentumMeans->cd(4); + hFittedPyStdDevs->Draw("E1"); // "E1" draws error bars + cFittedMomentumMeans->SetGrid(); + cFittedMomentumMeans->Update(); + // Save the canvas as a PNG file + cFittedMomentumMeans->SaveAs("fitted_momentum_means_stddevs.png"); + + + // ----------------------------------------------------------------------------- + // Create histograms of the beampipe parameters + // ----------------------------------------------------------------------------- + TH1F* hPipeRadii = CreateFittedHistogram("hPipeRadii", + "Radius [cm]", + pipeRadii, + {}, + "Pipe ID"); + TH1F* hPipeXPos = CreateFittedHistogram("hPipeXPos", + "X Position [cm]", + pipeXPos, + {}, + "Pipe ID"); + TH1F* hPipeZPos = CreateFittedHistogram("hPipeZPos", + "Z Position [cm]", + pipeZPos, + {}, + "Pipe ID"); + TH1F* hPipeRotations = CreateFittedHistogram("hPipeRotations", + "Rotation [rad]", + pipeRotation, + {}, + "Pipe ID"); + + // Create a canvas for the pipe parameters + TCanvas *cPipeParams = new TCanvas("cPipeParams", "Pipe Parameters", 1200, 400); + cPipeParams->Divide(4, 1); + cPipeParams->cd(1); + hPipeRadii->Draw(""); + cPipeParams->cd(2); + hPipeXPos->Draw(""); + cPipeParams->cd(3); + hPipeZPos->Draw(""); + cPipeParams->cd(4); + hPipeRotations->Draw(""); + cPipeParams->SetGrid(); + cPipeParams->Update(); + // Save the canvas as a PNG file + cPipeParams->SaveAs("pipe_parameters.png"); + + TFile *f = new TFile(outFile,"RECREATE"); + cXY->Write(); cX->Write(); cY->Write(); + cFittedMeans->Write(); + cFittedMomentumMeans->Write(); + cPipeParams->Write(); f->Close(); std::cout << "Saving events to file" << std::endl; - // ROOT::RDF::RSnapshotOptions opts; - // opts.fMode = "UPDATE"; - // d1.Snapshot("events",outFile,{"pipeID","endID","xpos","ypos","zpos","xmom","ymom","zmom","xpos_global","ypos_global","zpos_global","px_global","py_global","pz_global"},opts); + ROOT::RDF::RSnapshotOptions opts; + // opts.fMode = "UPDATE"; + // d1.Snapshot("events",outFile,{"pipeID","endID","pipeRadius","xpos","ypos","zpos","xmom","ymom","zmom","xpos_global","ypos_global","zpos_global","px_global","py_global","pz_global"},opts); } diff --git a/benchmarks/beamline/beamGPS.mac b/benchmarks/beamline/beamGPS.mac index 28d771ba..8ecd6bbb 100644 --- a/benchmarks/beamline/beamGPS.mac +++ b/benchmarks/beamline/beamGPS.mac @@ -4,7 +4,8 @@ /gps/pos/type Beam /gps/pos/sigma_x 0.119 mm /gps/pos/sigma_y 0.0107 mm -/gps/energy 18 GeV +/gps/energy 17.846 GeV +#/gps/energy 18 GeV /gps/particle e- /run/beamOn 100000 \ No newline at end of file diff --git a/benchmarks/beamline/functors.h b/benchmarks/beamline/functors.h index ee2414d0..fede4b84 100644 --- a/benchmarks/beamline/functors.h +++ b/benchmarks/beamline/functors.h @@ -100,4 +100,79 @@ struct globalToLocal{ private: dd4hep::Detector& detector; dd4hep::VolumeManager volumeManager; -}; \ No newline at end of file +}; + +struct volParams{ + double radius; + double xPos; + double yPos; + double zPos; + double rotation; +}; + +// Functor to get the volume element parameters from the CellID +struct getVolumeParametersFromCellID { + getVolumeParametersFromCellID(dd4hep::Detector& det) : detector(det) { + volumeManager = dd4hep::VolumeManager::getVolumeManager(det); + } + + ROOT::VecOps::RVec operator()(const RVecHits& evt) { + ROOT::VecOps::RVec result; + // Look up the detector element using the CellID + for(const auto& h : evt) { + auto cellID = h.cellID; + auto detelement = volumeManager.lookupDetElement(cellID); + const TGeoMatrix& transform = detelement.nominal().worldTransformation(); + const Double_t* translation = transform.GetTranslation(); + 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(); + volParams params{ + cone.rMax1(), + translation[0], + translation[1], + translation[2], + rotationAngleY + }; + result.push_back(params); + } + return result; + } + +private: + dd4hep::Detector& detector; + dd4hep::VolumeManager volumeManager; +}; + +TH1F* CreateFittedHistogram(const std::string& histName, + const std::string& histTitle, + const std::map& valueMap, + const std::map& errorMap, + const std::string& xAxisLabel) { + // Number of bins corresponds to the number of entries in the value map + int nPoints = valueMap.size(); + TH1F* hist = new TH1F(histName.c_str(), (";" + xAxisLabel + ";" + histTitle).c_str(), nPoints, 0.5, nPoints + 0.5); + + // Fill the histogram with values and errors, and set custom bin labels + int binIndex = 1; // Start from bin 1 + for (const auto& [name, value] : valueMap) { + hist->SetBinContent(binIndex, value); // Set the bin content + if(errorMap.size()==valueMap.size()){ + hist->SetBinError(binIndex, errorMap.at(name)); // Set the bin error + } + else{ + hist->SetBinError(binIndex, 0); // Set the bin error to 0 if not provided + } + hist->GetXaxis()->SetBinLabel(binIndex, name); // Set the bin label + binIndex++; + } + + // Customize the histogram + hist->SetMarkerStyle(20); + hist->SetMarkerSize(1.0); + hist->SetLineColor(kBlue); + hist->SetMarkerColor(kRed); + + return hist; +} \ No newline at end of file From 6329372609f67c3a5035d03686b71ed62a768c22 Mon Sep 17 00:00:00 2001 From: simonge Date: Fri, 2 May 2025 10:23:14 +0100 Subject: [PATCH 06/36] Add phasespace scripts --- benchmarks/beamline/Snakefile | 31 +- benchmarks/beamline/beamlineAnalysis.C | 453 ++++++++++++++++++ .../beamline/{beamGPS.mac => beamlineGPS.mac} | 0 .../{analysis.C => phasespaceAnalysis.C} | 2 +- benchmarks/beamline/phasespaceGPS.mac | 11 + 5 files changed, 489 insertions(+), 8 deletions(-) create mode 100644 benchmarks/beamline/beamlineAnalysis.C rename benchmarks/beamline/{beamGPS.mac => beamlineGPS.mac} (100%) rename benchmarks/beamline/{analysis.C => phasespaceAnalysis.C} (99%) create mode 100644 benchmarks/beamline/phasespaceGPS.mac diff --git a/benchmarks/beamline/Snakefile b/benchmarks/beamline/Snakefile index 13cd66a7..dc091e03 100644 --- a/benchmarks/beamline/Snakefile +++ b/benchmarks/beamline/Snakefile @@ -1,10 +1,10 @@ rule beamline_sim: input: - macro="beamGPS.mac", + macro="{test}GPS.mac", output: - "/scratch/EIC/G4out/beamline/beamlineTest{tag}.edm4hep.root", + "/scratch/EIC/G4out/beamline/{test}Test{tag}.edm4hep.root", log: - "/scratch/EIC/G4out/beamline/beamlineTest{tag}.edm4hep.root.log", + "/scratch/EIC/G4out/beamline/{test}Test{tag}.edm4hep.root.log", shell: """ exec npsim \ @@ -15,16 +15,33 @@ rule beamline_sim: --outputFile {output} \ """ +rule beamline_geometry_image: + output: + "geometry_image{tag}.png" + params: + xml=os.getenv("DETECTOR_PATH")+"/epic_ip6_extended.xml", + shell: + """ + exec npsim \ + --compactFile {params.xml} \ + --runType vis \ + --enableG4GPS \ + --macroFile visbeam.mac + mv geometry_image.png {output} + """ + rule beamline_analysis: input: - "/scratch/EIC/G4out/beamline/beamlineTest{tag}.edm4hep.root", + data="/scratch/EIC/G4out/beamline/{test}Test{tag}.edm4hep.root", + script="{test}Analysis.C", output: - "/scratch/EIC/ReconOut/beamline/beamlineTestAnalysis{tag}.root" + "/scratch/EIC/ReconOut/beamline/{test}TestAnalysis{tag}.root" log: - "/scratch/EIC/ReconOut/beamline/beamlineTestAnalysis{tag}.root.log" + "/scratch/EIC/ReconOut/beamline/{test}TestAnalysis{tag}.root.log" params: xml=os.getenv("DETECTOR_PATH")+"/epic_ip6_extended.xml", shell: """ - root -l -b -q 'analysis.C("{input}", "{output}", "{params.xml}")' + root -l -b -q 'analysis.C("{input.data}", "{output}", "{params.xml}")' """ + diff --git a/benchmarks/beamline/beamlineAnalysis.C b/benchmarks/beamline/beamlineAnalysis.C new file mode 100644 index 00000000..abd677a8 --- /dev/null +++ b/benchmarks/beamline/beamlineAnalysis.C @@ -0,0 +1,453 @@ +// Script to plot the x and y positions and momentum of beamline particles as they pass through the magnets + +#include +#include "DD4hep/Detector.h" +#include "DDRec/CellIDPositionConverter.h" +#include "edm4hep/MCParticleCollection.h" +#include "edm4hep/SimTrackerHitCollection.h" +#include "edm4hep/SimCalorimeterHitCollection.h" +#include "edm4hep/Vector3f.h" +#include "edm4hep/Vector3d.h" +#include "ROOT/RDataFrame.hxx" +#include "ROOT/RDF/RInterface.hxx" +#include "ROOT/RVec.hxx" +#include "functors.h" +#include "TCanvas.h" +#include "TStyle.h" + + +using RVecS = ROOT::VecOps::RVec; +using RNode = ROOT::RDF::RNode; + +void beamlineAnalysis( TString inFile = "/scratch/EIC/G4out/beamline/beamlineTest.edm4hep.root", + TString outFile = "output.root", + std::string compactName = "/home/simong/EIC/epic/install/share/epic/epic_ip6_extended.xml"){ + + //Set ROOT style + gStyle->SetPadLeftMargin(0.1); // Set left margin + gStyle->SetPadRightMargin(0.0); // Set right margin + gStyle->SetPadTopMargin(0.0); // Set top margin + gStyle->SetPadBottomMargin(0.1);// Set bottom margin + gStyle->SetTitleAlign(13); + gStyle->SetTitleX(0.12); // Place the title on the top right + gStyle->SetTitleY(0.985); // Place the title on the top right + gStyle->SetTitleSize(0.08, "t"); + gStyle->SetTitleXOffset(1.0); // Adjust y-axis title offset + gStyle->SetTitleYOffset(1.0); // Adjust y-axis title offset + gStyle->SetOptStat(0); + + //Set implicit multi-threading + ROOT::EnableImplicitMT(); + + //Load the detector config + dd4hep::Detector& detector = dd4hep::Detector::getInstance(); + detector.fromCompact(compactName); + + ROOT::RDataFrame d0("events",inFile, {"BackwardsBeamlineHits"}); + RNode d1 = d0; + RVecS colNames = d1.GetColumnNames(); + + //Set number of entries to process + // d1 = d1.Range(0,1000); + + //Get the collection + std::string readoutName = "BackwardsBeamlineHits"; + + std::cout << "Running lazy RDataframe execution" << std::endl; + + if(Any(colNames==readoutName)){ + + d1 = d1.Define("pipeID",getSubID("pipe",detector),{readoutName}) + .Define("endID",getSubID("end",detector),{readoutName}) + .Define("pipeParameters",getVolumeParametersFromCellID(detector),{readoutName}) + .Define("pipeRadius",[](const ROOT::VecOps::RVec& params) { + ROOT::VecOps::RVec radii; + for (const auto& param : params) { + radii.push_back(param.radius); + } + return radii; + }, {"pipeParameters"}) + .Define("xdet",[](const ROOT::VecOps::RVec& params) { + ROOT::VecOps::RVec xPos; + for (const auto& param : params) { + xPos.push_back(param.xPos); + } + return xPos; + }, {"pipeParameters"}) + .Define("zdet",[](const ROOT::VecOps::RVec& params) { + ROOT::VecOps::RVec zPos; + for (const auto& param : params) { + zPos.push_back(param.zPos); + } + return zPos; + }, {"pipeParameters"}) + .Define("rotation",[](const ROOT::VecOps::RVec& params) { + ROOT::VecOps::RVec rotation; + for (const auto& param : params) { + rotation.push_back(param.rotation); + } + return rotation; + }, {"pipeParameters"}); + + + //global x,y,z position and momentum + d1 = d1.Define("xpos_global","BackwardsBeamlineHits.position.x") + .Define("ypos_global","BackwardsBeamlineHits.position.y") + .Define("zpos_global","BackwardsBeamlineHits.position.z") + .Define("px_global","BackwardsBeamlineHits.momentum.x") + .Define("py_global","BackwardsBeamlineHits.momentum.y") + .Define("pz_global","BackwardsBeamlineHits.momentum.z"); + + d1 = d1.Define("hitPosMom",globalToLocal(detector),{readoutName}) + .Define("xpos","hitPosMom[0]") + .Define("ypos","hitPosMom[1]") + .Define("zpos","hitPosMom[2]") + .Define("xmomMag","hitPosMom[3]") + .Define("ymomMag","hitPosMom[4]") + .Define("zmomMag","hitPosMom[5]") + .Define("momMag","sqrt(xmomMag*xmomMag+ymomMag*ymomMag+zmomMag*zmomMag)") + .Define("xmom","xmomMag/momMag") + .Define("ymom","ymomMag/momMag") + .Define("zmom","zmomMag/momMag"); + + } + else{ + std::cout << "Collection " << readoutName << " not found in file" << std::endl; + return; + } + + // Calculate the maximum pipe radius for plotting + auto maxPipeRadius = 1.2*d1.Max("pipeRadius").GetValue(); + + std::cout << "Executing Analysis and creating histograms" << std::endl; + + //Create array of histogram results + std::map> hHistsxy; + std::map> hHistsxyZoom; + std::map> hHistsxpx; + std::map> hHistsypy; + + std::map xMeans; + std::map yMeans; + std::map xStdDevs; + std::map yStdDevs; + std::map pxMeans; + std::map pyMeans; + std::map pxStdDevs; + std::map pyStdDevs; + + //Fit paremeter and error maps + std::map xMeanFit; + std::map yMeanFit; + std::map xMeanFitErr; + std::map yMeanFitErr; + std::map xStdDevFit; + std::map yStdDevFit; + std::map xStdDevFitErr; + std::map yStdDevFitErr; + std::map pxMeanFit; + std::map pyMeanFit; + std::map pxMeanFitErr; + std::map pyMeanFitErr; + std::map pxStdDevFit; + std::map pyStdDevFit; + std::map pxStdDevFitErr; + std::map pyStdDevFitErr; + + std::map pipeRadii; + std::map pipeXPos; + std::map pipeZPos; + std::map pipeRotation; + + auto xmin = d1.Min("xpos").GetValue(); + auto xmax = d1.Max("xpos").GetValue(); + auto ymin = d1.Min("ypos").GetValue(); + auto ymax = d1.Max("ypos").GetValue(); + auto pxmin = d1.Min("xmom").GetValue(); + auto pxmax = d1.Max("xmom").GetValue(); + auto pymin = d1.Min("ymom").GetValue(); + auto pymax = d1.Max("ymom").GetValue(); + + //Create histograms + for(int i=0; i<=7; i++){ + + std::string name = "pipeID"; + name += std::to_string(i); + auto filterDF = d1.Define("xposf","xpos["+std::to_string(i)+"]") + .Define("yposf","ypos["+std::to_string(i)+"]") + .Define("xmomf","xmom["+std::to_string(i)+"]") + .Define("ymomf","ymom["+std::to_string(i)+"]") + .Define("pipeRadiusf","pipeRadius["+std::to_string(i)+"]") + .Define("xdetf","xdet["+std::to_string(i)+"]") + .Define("zdetf","zdet["+std::to_string(i)+"]") + .Define("rotationf","rotation["+std::to_string(i)+"]"); + + + //Calculate Min and Max values + auto xminf = filterDF.Min("xposf").GetValue(); + auto xmaxf = filterDF.Max("xposf").GetValue(); + auto yminf = filterDF.Min("yposf").GetValue(); + auto ymaxf = filterDF.Max("yposf").GetValue(); + auto pxminf = filterDF.Min("xmomf").GetValue(); + auto pxmaxf = filterDF.Max("xmomf").GetValue(); + auto pyminf = filterDF.Min("ymomf").GetValue(); + auto pymaxf = filterDF.Max("ymomf").GetValue(); + // Calculate means and standard deviations + xMeans[name] = filterDF.Mean("xposf").GetValue(); + yMeans[name] = filterDF.Mean("yposf").GetValue(); + xStdDevs[name] = filterDF.StdDev("xposf").GetValue(); + yStdDevs[name] = filterDF.StdDev("yposf").GetValue(); + pxMeans[name] = filterDF.Mean("xmomf").GetValue(); + pyMeans[name] = filterDF.Mean("ymomf").GetValue(); + pxStdDevs[name] = filterDF.StdDev("xmomf").GetValue(); + pyStdDevs[name] = filterDF.StdDev("ymomf").GetValue(); + + // Calculate axes for zoomed beamspot histogram so that it is quare around the mean x and y + double halfrange = std::max({xMeans[name]-xminf, xmaxf-xMeans[name], yMeans[name]-yminf, ymaxf-yMeans[name]}); + double xMinZoom = xMeans[name] - halfrange; + double xMaxZoom = xMeans[name] + halfrange; + double yMinZoom = yMeans[name] - halfrange; + double yMaxZoom = yMeans[name] + halfrange; + + TString beamspotName = "Beamspot ID"+std::to_string(i)+";x offset [cm]; y offset [cm]"; + TString xyname = name+";x Offset [cm]; y Offset [cm]"; + TString xname = name+";x Offset [cm]; x trajectory component"; + TString yname = name+";y Offset [cm]; y trajectory component"; + hHistsxy[name] = filterDF.Histo2D({beamspotName,beamspotName,400,-maxPipeRadius,maxPipeRadius,400,-maxPipeRadius,maxPipeRadius},"xposf","yposf"); + hHistsxyZoom[name] = filterDF.Histo2D({xyname,xyname,100,xMinZoom,xMaxZoom,100,yMinZoom,yMaxZoom},"xposf","yposf"); + hHistsxpx[name] = filterDF.Histo2D({xname,xname,400,xmin,xmax,400,pxmin,pxmax},"xposf","xmomf"); + hHistsypy[name] = filterDF.Histo2D({yname,yname,400,ymin,ymax,400,pymin,pymax},"yposf","ymomf"); + + //Parameters of the pipe + pipeRadii[name] = filterDF.Max("pipeRadiusf").GetValue(); + pipeXPos[name] = filterDF.Max("xdetf").GetValue(); + pipeZPos[name] = filterDF.Max("zdetf").GetValue(); + pipeRotation[name] = filterDF.Max("rotationf").GetValue(); + + //Fit gaussian to the x, y, px and py distributions + auto xhist = hHistsxy[name]->ProjectionX(); + auto yhist = hHistsxy[name]->ProjectionY(); + auto pxhist = hHistsxpx[name]->ProjectionY(); + auto pyhist = hHistsypy[name]->ProjectionY(); + xhist->Fit("gaus","Q"); + yhist->Fit("gaus","Q"); + pxhist->Fit("gaus","Q"); + pyhist->Fit("gaus","Q"); + //Get the fit parameters and errors + xMeanFit[name] = xhist->GetFunction("gaus")->GetParameter(1); + yMeanFit[name] = yhist->GetFunction("gaus")->GetParameter(1); + xMeanFitErr[name] = xhist->GetFunction("gaus")->GetParError(1); + yMeanFitErr[name] = yhist->GetFunction("gaus")->GetParError(1); + xStdDevFit[name] = xhist->GetFunction("gaus")->GetParameter(2); + yStdDevFit[name] = yhist->GetFunction("gaus")->GetParameter(2); + xStdDevFitErr[name] = xhist->GetFunction("gaus")->GetParError(2); + yStdDevFitErr[name] = yhist->GetFunction("gaus")->GetParError(2); + pxMeanFit[name] = pxhist->GetFunction("gaus")->GetParameter(1); + pyMeanFit[name] = pyhist->GetFunction("gaus")->GetParameter(1); + pxMeanFitErr[name] = pxhist->GetFunction("gaus")->GetParError(1); + pyMeanFitErr[name] = pyhist->GetFunction("gaus")->GetParError(1); + pxStdDevFit[name] = pxhist->GetFunction("gaus")->GetParameter(2); + pyStdDevFit[name] = pyhist->GetFunction("gaus")->GetParameter(2); + pxStdDevFitErr[name] = pxhist->GetFunction("gaus")->GetParError(2); + pyStdDevFitErr[name] = pyhist->GetFunction("gaus")->GetParError(2); + + } + + // Create histograms of the beamspot + TCanvas *cXY = new TCanvas("beamspot_canvas","beamspot_canvas",3000,1600); + cXY->Divide(4,2); + int i=1; + for(auto [name,h] : hHistsxy){ + // Get the pipe radius for this histogram + auto pipeRadius = pipeRadii[name]; + 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"); + + // Add zoomed version in the top-right corner + TPad *pad = new TPad("zoomPad", "Zoomed View", 0.65, 0.65, 1.0, 1.0); + pad->SetFillStyle(0); // Transparent background + pad->Draw(); + pad->cd(); + // Draw the zoomed histogram without its title or axis titles + hHistsxyZoom[name]->SetTitle(""); + hHistsxyZoom[name]->GetXaxis()->SetTitle(""); + hHistsxyZoom[name]->GetYaxis()->SetTitle(""); + hHistsxyZoom[name]->Draw("col"); + cXY->cd(); // Return to the main canvas + } + + // x-px canvas + TCanvas *cX = new TCanvas("x_px_canvas","x_px_canvas",3000,1600); + cX->Divide(4,2); + + i=1; + //Write histograms to file + for(auto& h : hHistsxpx){ + cX->cd(i++); + h.second->Draw("col"); + } + + // y-py canvas + TCanvas *cY = new TCanvas("y_py_canvas","y_py_canvas",3000,1600); + cY->Divide(4,2); + + i=1; + for(auto& h : hHistsypy){ + cY->cd(i++); + h.second->Draw("col"); + } + + // Save 2D canvases as pngs + cXY->SaveAs("beamspot.png"); + cX->SaveAs("x_px.png"); + cY->SaveAs("y_py.png"); + + // --------------------------------------------------------------------------- + // Create histograms showing the fitted means and standard deviations of the positions and momenta + // --------------------------------------------------------------------------- + + // Create histograms for fitted X means and standard deviations + TH1F* hFittedXMeans = CreateFittedHistogram("hFittedXMeans", + "Mean X Offset [cm]", + xMeanFit, + xMeanFitErr, + "Pipe ID"); + + TH1F* hFittedXStdDevs = CreateFittedHistogram("hFittedXStdDevs", + "Std Deviation X Offset [cm]", + xStdDevFit, + xStdDevFitErr, + "Pipe ID"); + + // Create histograms for fitted Y means and standard deviations + TH1F* hFittedYMeans = CreateFittedHistogram("hFittedYMeans", + "Mean Y Offset [cm]", + yMeanFit, + yMeanFitErr, + "Pipe ID"); + + TH1F* hFittedYStdDevs = CreateFittedHistogram("hFittedYStdDevs", + "Std Deviation Y Offset [cm]", + yStdDevFit, + yStdDevFitErr, + "Pipe ID"); + + TH1F* hFittedPxMeans = CreateFittedHistogram("hFittedPxMeans", + "Mean Px", + pxMeanFit, + pxMeanFitErr, + "Pipe ID"); + TH1F* hFittedPyMeans = CreateFittedHistogram("hFittedPyMeans", + "Mean Py", + pyMeanFit, + pyMeanFitErr, + "Pipe ID"); + TH1F* hFittedPxStdDevs = CreateFittedHistogram("hFittedPxStdDevs", + "Std Deviation Px", + pxStdDevFit, + pxStdDevFitErr, + "Pipe ID"); + TH1F* hFittedPyStdDevs = CreateFittedHistogram("hFittedPyStdDevs", + "Std Deviation Py", + pyStdDevFit, + pyStdDevFitErr, + "Pipe ID"); + + + // Create a canvas for the fitted beamspot means and standard deviations + TCanvas *cFittedMeans = new TCanvas("cFittedMeans", "Fitted Beamspot Means and Std Deviation", 1200, 800); + cFittedMeans->Divide(2, 2); + cFittedMeans->cd(1); + hFittedXMeans->Draw("E1"); // "E1" draws error bars + cFittedMeans->cd(2); + hFittedXStdDevs->Draw("E1"); // "E1" draws error bars + cFittedMeans->cd(3); + hFittedYMeans->Draw("E1"); // "E1" draws error bars + cFittedMeans->cd(4); + hFittedYStdDevs->Draw("E1"); // "E1" draws error bars + cFittedMeans->SetGrid(); + cFittedMeans->Update(); + // Save the canvas as a PNG file + cFittedMeans->SaveAs("fitted_means_stddevs.png"); + + // Create a canvas for the fitted momentum means and standard deviations + TCanvas *cFittedMomentumMeans = new TCanvas("cFittedMomentumMeans", "Fitted Momentum Means and Std Deviation", 1200, 800); + cFittedMomentumMeans->Divide(2, 2); + cFittedMomentumMeans->cd(1); + hFittedPxMeans->Draw("E1"); // "E1" draws error bars + cFittedMomentumMeans->cd(2); + hFittedPxStdDevs->Draw("E1"); // "E1" draws error bars + cFittedMomentumMeans->cd(3); + hFittedPyMeans->Draw("E1"); // "E1" draws error bars + cFittedMomentumMeans->cd(4); + hFittedPyStdDevs->Draw("E1"); // "E1" draws error bars + cFittedMomentumMeans->SetGrid(); + cFittedMomentumMeans->Update(); + // Save the canvas as a PNG file + cFittedMomentumMeans->SaveAs("fitted_momentum_means_stddevs.png"); + + + // ----------------------------------------------------------------------------- + // Create histograms of the beampipe parameters + // ----------------------------------------------------------------------------- + TH1F* hPipeRadii = CreateFittedHistogram("hPipeRadii", + "Radius [cm]", + pipeRadii, + {}, + "Pipe ID"); + TH1F* hPipeXPos = CreateFittedHistogram("hPipeXPos", + "X Position [cm]", + pipeXPos, + {}, + "Pipe ID"); + TH1F* hPipeZPos = CreateFittedHistogram("hPipeZPos", + "Z Position [cm]", + pipeZPos, + {}, + "Pipe ID"); + TH1F* hPipeRotations = CreateFittedHistogram("hPipeRotations", + "Rotation [rad]", + pipeRotation, + {}, + "Pipe ID"); + + // Create a canvas for the pipe parameters + TCanvas *cPipeParams = new TCanvas("cPipeParams", "Pipe Parameters", 1200, 400); + cPipeParams->Divide(4, 1); + cPipeParams->cd(1); + hPipeRadii->Draw(""); + cPipeParams->cd(2); + hPipeXPos->Draw(""); + cPipeParams->cd(3); + hPipeZPos->Draw(""); + cPipeParams->cd(4); + hPipeRotations->Draw(""); + cPipeParams->SetGrid(); + cPipeParams->Update(); + // Save the canvas as a PNG file + cPipeParams->SaveAs("pipe_parameters.png"); + + + TFile *f = new TFile(outFile,"RECREATE"); + cXY->Write(); + cX->Write(); + cY->Write(); + cFittedMeans->Write(); + cFittedMomentumMeans->Write(); + cPipeParams->Write(); + + f->Close(); + + std::cout << "Saving events to file" << std::endl; + + ROOT::RDF::RSnapshotOptions opts; + // opts.fMode = "UPDATE"; + // d1.Snapshot("events",outFile,{"pipeID","endID","pipeRadius","xpos","ypos","zpos","xmom","ymom","zmom","xpos_global","ypos_global","zpos_global","px_global","py_global","pz_global"},opts); +} diff --git a/benchmarks/beamline/beamGPS.mac b/benchmarks/beamline/beamlineGPS.mac similarity index 100% rename from benchmarks/beamline/beamGPS.mac rename to benchmarks/beamline/beamlineGPS.mac diff --git a/benchmarks/beamline/analysis.C b/benchmarks/beamline/phasespaceAnalysis.C similarity index 99% rename from benchmarks/beamline/analysis.C rename to benchmarks/beamline/phasespaceAnalysis.C index f37c2ce1..770e52cd 100644 --- a/benchmarks/beamline/analysis.C +++ b/benchmarks/beamline/phasespaceAnalysis.C @@ -19,7 +19,7 @@ using RVecS = ROOT::VecOps::RVec; using RNode = ROOT::RDF::RNode; -void analysis( TString inFile = "/scratch/EIC/G4out/beamline/beamlineTest.edm4hep.root", +void phasespaceAnalysis( TString inFile = "/scratch/EIC/G4out/beamline/beamlineTest.edm4hep.root", TString outFile = "output.root", std::string compactName = "/home/simong/EIC/epic/install/share/epic/epic_ip6_extended.xml"){ diff --git a/benchmarks/beamline/phasespaceGPS.mac b/benchmarks/beamline/phasespaceGPS.mac new file mode 100644 index 00000000..9ef8ed23 --- /dev/null +++ b/benchmarks/beamline/phasespaceGPS.mac @@ -0,0 +1,11 @@ +/gps/ang/type beam2d +/gps/ang/sigma_x 0.011 rad # 11 mrad +/gps/ang/sigma_y 0.011 rad # 11 mrad +/gps/pos/type Beam +/gps/pos/sigma_x 0.119 mm +/gps/pos/sigma_y 0.0107 mm +/gps/energy 17.846 GeV +#/gps/energy 18 GeV +/gps/particle e- + +/run/beamOn 100000 \ No newline at end of file From 67d765c8a33a9c808e75e41f9411bed67bae35d9 Mon Sep 17 00:00:00 2001 From: simonge Date: Fri, 2 May 2025 14:16:13 +0100 Subject: [PATCH 07/36] Change name of functors --- benchmarks/beamline/beamlineAnalysis.C | 2 +- benchmarks/beamline/phasespaceAnalysis.C | 14 +++++++------- benchmarks/beamline/phasespaceGPS.mac | 2 +- .../beamline/{functors.h => shared_functions.h} | 0 4 files changed, 9 insertions(+), 9 deletions(-) rename benchmarks/beamline/{functors.h => shared_functions.h} (100%) diff --git a/benchmarks/beamline/beamlineAnalysis.C b/benchmarks/beamline/beamlineAnalysis.C index abd677a8..32cd51cd 100644 --- a/benchmarks/beamline/beamlineAnalysis.C +++ b/benchmarks/beamline/beamlineAnalysis.C @@ -11,7 +11,7 @@ #include "ROOT/RDataFrame.hxx" #include "ROOT/RDF/RInterface.hxx" #include "ROOT/RVec.hxx" -#include "functors.h" +#include "shared_functions.h" #include "TCanvas.h" #include "TStyle.h" diff --git a/benchmarks/beamline/phasespaceAnalysis.C b/benchmarks/beamline/phasespaceAnalysis.C index 770e52cd..6ee5b05c 100644 --- a/benchmarks/beamline/phasespaceAnalysis.C +++ b/benchmarks/beamline/phasespaceAnalysis.C @@ -11,7 +11,7 @@ #include "ROOT/RDataFrame.hxx" #include "ROOT/RDF/RInterface.hxx" #include "ROOT/RVec.hxx" -#include "functors.h" +#include "shared_functions.h" #include "TCanvas.h" #include "TStyle.h" @@ -305,9 +305,9 @@ void phasespaceAnalysis( TString inFile = "/scratch/EIC/G4out/beamline/bea } // Save 2D canvases as pngs - cXY->SaveAs("beamspot.png"); - cX->SaveAs("x_px.png"); - cY->SaveAs("y_py.png"); + // cXY->SaveAs("beamspot.png"); + // cX->SaveAs("x_px.png"); + // cY->SaveAs("y_py.png"); // --------------------------------------------------------------------------- // Create histograms showing the fitted means and standard deviations of the positions and momenta @@ -375,7 +375,7 @@ void phasespaceAnalysis( TString inFile = "/scratch/EIC/G4out/beamline/bea cFittedMeans->SetGrid(); cFittedMeans->Update(); // Save the canvas as a PNG file - cFittedMeans->SaveAs("fitted_means_stddevs.png"); + // cFittedMeans->SaveAs("fitted_means_stddevs.png"); // Create a canvas for the fitted momentum means and standard deviations TCanvas *cFittedMomentumMeans = new TCanvas("cFittedMomentumMeans", "Fitted Momentum Means and Std Deviation", 1200, 800); @@ -391,7 +391,7 @@ void phasespaceAnalysis( TString inFile = "/scratch/EIC/G4out/beamline/bea cFittedMomentumMeans->SetGrid(); cFittedMomentumMeans->Update(); // Save the canvas as a PNG file - cFittedMomentumMeans->SaveAs("fitted_momentum_means_stddevs.png"); + // cFittedMomentumMeans->SaveAs("fitted_momentum_means_stddevs.png"); // ----------------------------------------------------------------------------- @@ -432,7 +432,7 @@ void phasespaceAnalysis( TString inFile = "/scratch/EIC/G4out/beamline/bea cPipeParams->SetGrid(); cPipeParams->Update(); // Save the canvas as a PNG file - cPipeParams->SaveAs("pipe_parameters.png"); + // cPipeParams->SaveAs("pipe_parameters.png"); TFile *f = new TFile(outFile,"RECREATE"); diff --git a/benchmarks/beamline/phasespaceGPS.mac b/benchmarks/beamline/phasespaceGPS.mac index 9ef8ed23..5d33ceff 100644 --- a/benchmarks/beamline/phasespaceGPS.mac +++ b/benchmarks/beamline/phasespaceGPS.mac @@ -8,4 +8,4 @@ #/gps/energy 18 GeV /gps/particle e- -/run/beamOn 100000 \ No newline at end of file +/run/beamOn 1000000 \ No newline at end of file diff --git a/benchmarks/beamline/functors.h b/benchmarks/beamline/shared_functions.h similarity index 100% rename from benchmarks/beamline/functors.h rename to benchmarks/beamline/shared_functions.h From 3afabcb8d8f3e470c958d052a716537b483d7049 Mon Sep 17 00:00:00 2001 From: simonge Date: Wed, 4 Jun 2025 14:53:53 +0100 Subject: [PATCH 08/36] Remove image and acceptance benchmarks for now --- benchmarks/beamline/Snakefile | 15 - benchmarks/beamline/phasespaceAnalysis.C | 453 ----------------------- benchmarks/beamline/phasespaceGPS.mac | 11 - 3 files changed, 479 deletions(-) delete mode 100644 benchmarks/beamline/phasespaceAnalysis.C delete mode 100644 benchmarks/beamline/phasespaceGPS.mac diff --git a/benchmarks/beamline/Snakefile b/benchmarks/beamline/Snakefile index dc091e03..88de43c1 100644 --- a/benchmarks/beamline/Snakefile +++ b/benchmarks/beamline/Snakefile @@ -15,21 +15,6 @@ rule beamline_sim: --outputFile {output} \ """ -rule beamline_geometry_image: - output: - "geometry_image{tag}.png" - params: - xml=os.getenv("DETECTOR_PATH")+"/epic_ip6_extended.xml", - shell: - """ - exec npsim \ - --compactFile {params.xml} \ - --runType vis \ - --enableG4GPS \ - --macroFile visbeam.mac - mv geometry_image.png {output} - """ - rule beamline_analysis: input: data="/scratch/EIC/G4out/beamline/{test}Test{tag}.edm4hep.root", diff --git a/benchmarks/beamline/phasespaceAnalysis.C b/benchmarks/beamline/phasespaceAnalysis.C deleted file mode 100644 index 6ee5b05c..00000000 --- a/benchmarks/beamline/phasespaceAnalysis.C +++ /dev/null @@ -1,453 +0,0 @@ -// Script to plot the x and y positions and momentum of beamline particles as they pass through the magnets - -#include -#include "DD4hep/Detector.h" -#include "DDRec/CellIDPositionConverter.h" -#include "edm4hep/MCParticleCollection.h" -#include "edm4hep/SimTrackerHitCollection.h" -#include "edm4hep/SimCalorimeterHitCollection.h" -#include "edm4hep/Vector3f.h" -#include "edm4hep/Vector3d.h" -#include "ROOT/RDataFrame.hxx" -#include "ROOT/RDF/RInterface.hxx" -#include "ROOT/RVec.hxx" -#include "shared_functions.h" -#include "TCanvas.h" -#include "TStyle.h" - - -using RVecS = ROOT::VecOps::RVec; -using RNode = ROOT::RDF::RNode; - -void phasespaceAnalysis( TString inFile = "/scratch/EIC/G4out/beamline/beamlineTest.edm4hep.root", - TString outFile = "output.root", - std::string compactName = "/home/simong/EIC/epic/install/share/epic/epic_ip6_extended.xml"){ - - //Set ROOT style - gStyle->SetPadLeftMargin(0.1); // Set left margin - gStyle->SetPadRightMargin(0.0); // Set right margin - gStyle->SetPadTopMargin(0.0); // Set top margin - gStyle->SetPadBottomMargin(0.1);// Set bottom margin - gStyle->SetTitleAlign(13); - gStyle->SetTitleX(0.12); // Place the title on the top right - gStyle->SetTitleY(0.985); // Place the title on the top right - gStyle->SetTitleSize(0.08, "t"); - gStyle->SetTitleXOffset(1.0); // Adjust y-axis title offset - gStyle->SetTitleYOffset(1.0); // Adjust y-axis title offset - gStyle->SetOptStat(0); - - //Set implicit multi-threading - ROOT::EnableImplicitMT(); - - //Load the detector config - dd4hep::Detector& detector = dd4hep::Detector::getInstance(); - detector.fromCompact(compactName); - - ROOT::RDataFrame d0("events",inFile, {"BackwardsBeamlineHits"}); - RNode d1 = d0; - RVecS colNames = d1.GetColumnNames(); - - //Set number of entries to process - // d1 = d1.Range(0,1000); - - //Get the collection - std::string readoutName = "BackwardsBeamlineHits"; - - std::cout << "Running lazy RDataframe execution" << std::endl; - - if(Any(colNames==readoutName)){ - - d1 = d1.Define("pipeID",getSubID("pipe",detector),{readoutName}) - .Define("endID",getSubID("end",detector),{readoutName}) - .Define("pipeParameters",getVolumeParametersFromCellID(detector),{readoutName}) - .Define("pipeRadius",[](const ROOT::VecOps::RVec& params) { - ROOT::VecOps::RVec radii; - for (const auto& param : params) { - radii.push_back(param.radius); - } - return radii; - }, {"pipeParameters"}) - .Define("xdet",[](const ROOT::VecOps::RVec& params) { - ROOT::VecOps::RVec xPos; - for (const auto& param : params) { - xPos.push_back(param.xPos); - } - return xPos; - }, {"pipeParameters"}) - .Define("zdet",[](const ROOT::VecOps::RVec& params) { - ROOT::VecOps::RVec zPos; - for (const auto& param : params) { - zPos.push_back(param.zPos); - } - return zPos; - }, {"pipeParameters"}) - .Define("rotation",[](const ROOT::VecOps::RVec& params) { - ROOT::VecOps::RVec rotation; - for (const auto& param : params) { - rotation.push_back(param.rotation); - } - return rotation; - }, {"pipeParameters"}); - - - //global x,y,z position and momentum - d1 = d1.Define("xpos_global","BackwardsBeamlineHits.position.x") - .Define("ypos_global","BackwardsBeamlineHits.position.y") - .Define("zpos_global","BackwardsBeamlineHits.position.z") - .Define("px_global","BackwardsBeamlineHits.momentum.x") - .Define("py_global","BackwardsBeamlineHits.momentum.y") - .Define("pz_global","BackwardsBeamlineHits.momentum.z"); - - d1 = d1.Define("hitPosMom",globalToLocal(detector),{readoutName}) - .Define("xpos","hitPosMom[0]") - .Define("ypos","hitPosMom[1]") - .Define("zpos","hitPosMom[2]") - .Define("xmomMag","hitPosMom[3]") - .Define("ymomMag","hitPosMom[4]") - .Define("zmomMag","hitPosMom[5]") - .Define("momMag","sqrt(xmomMag*xmomMag+ymomMag*ymomMag+zmomMag*zmomMag)") - .Define("xmom","xmomMag/momMag") - .Define("ymom","ymomMag/momMag") - .Define("zmom","zmomMag/momMag"); - - } - else{ - std::cout << "Collection " << readoutName << " not found in file" << std::endl; - return; - } - - // Calculate the maximum pipe radius for plotting - auto maxPipeRadius = 1.2*d1.Max("pipeRadius").GetValue(); - - std::cout << "Executing Analysis and creating histograms" << std::endl; - - //Create array of histogram results - std::map> hHistsxy; - std::map> hHistsxyZoom; - std::map> hHistsxpx; - std::map> hHistsypy; - - std::map xMeans; - std::map yMeans; - std::map xStdDevs; - std::map yStdDevs; - std::map pxMeans; - std::map pyMeans; - std::map pxStdDevs; - std::map pyStdDevs; - - //Fit paremeter and error maps - std::map xMeanFit; - std::map yMeanFit; - std::map xMeanFitErr; - std::map yMeanFitErr; - std::map xStdDevFit; - std::map yStdDevFit; - std::map xStdDevFitErr; - std::map yStdDevFitErr; - std::map pxMeanFit; - std::map pyMeanFit; - std::map pxMeanFitErr; - std::map pyMeanFitErr; - std::map pxStdDevFit; - std::map pyStdDevFit; - std::map pxStdDevFitErr; - std::map pyStdDevFitErr; - - std::map pipeRadii; - std::map pipeXPos; - std::map pipeZPos; - std::map pipeRotation; - - auto xmin = d1.Min("xpos").GetValue(); - auto xmax = d1.Max("xpos").GetValue(); - auto ymin = d1.Min("ypos").GetValue(); - auto ymax = d1.Max("ypos").GetValue(); - auto pxmin = d1.Min("xmom").GetValue(); - auto pxmax = d1.Max("xmom").GetValue(); - auto pymin = d1.Min("ymom").GetValue(); - auto pymax = d1.Max("ymom").GetValue(); - - //Create histograms - for(int i=0; i<=7; i++){ - - std::string name = "pipeID"; - name += std::to_string(i); - auto filterDF = d1.Define("xposf","xpos["+std::to_string(i)+"]") - .Define("yposf","ypos["+std::to_string(i)+"]") - .Define("xmomf","xmom["+std::to_string(i)+"]") - .Define("ymomf","ymom["+std::to_string(i)+"]") - .Define("pipeRadiusf","pipeRadius["+std::to_string(i)+"]") - .Define("xdetf","xdet["+std::to_string(i)+"]") - .Define("zdetf","zdet["+std::to_string(i)+"]") - .Define("rotationf","rotation["+std::to_string(i)+"]"); - - - //Calculate Min and Max values - auto xminf = filterDF.Min("xposf").GetValue(); - auto xmaxf = filterDF.Max("xposf").GetValue(); - auto yminf = filterDF.Min("yposf").GetValue(); - auto ymaxf = filterDF.Max("yposf").GetValue(); - auto pxminf = filterDF.Min("xmomf").GetValue(); - auto pxmaxf = filterDF.Max("xmomf").GetValue(); - auto pyminf = filterDF.Min("ymomf").GetValue(); - auto pymaxf = filterDF.Max("ymomf").GetValue(); - // Calculate means and standard deviations - xMeans[name] = filterDF.Mean("xposf").GetValue(); - yMeans[name] = filterDF.Mean("yposf").GetValue(); - xStdDevs[name] = filterDF.StdDev("xposf").GetValue(); - yStdDevs[name] = filterDF.StdDev("yposf").GetValue(); - pxMeans[name] = filterDF.Mean("xmomf").GetValue(); - pyMeans[name] = filterDF.Mean("ymomf").GetValue(); - pxStdDevs[name] = filterDF.StdDev("xmomf").GetValue(); - pyStdDevs[name] = filterDF.StdDev("ymomf").GetValue(); - - // Calculate axes for zoomed beamspot histogram so that it is quare around the mean x and y - double halfrange = std::max({xMeans[name]-xminf, xmaxf-xMeans[name], yMeans[name]-yminf, ymaxf-yMeans[name]}); - double xMinZoom = xMeans[name] - halfrange; - double xMaxZoom = xMeans[name] + halfrange; - double yMinZoom = yMeans[name] - halfrange; - double yMaxZoom = yMeans[name] + halfrange; - - TString beamspotName = "Beamspot ID"+std::to_string(i)+";x offset [cm]; y offset [cm]"; - TString xyname = name+";x Offset [cm]; y Offset [cm]"; - TString xname = name+";x Offset [cm]; x trajectory component"; - TString yname = name+";y Offset [cm]; y trajectory component"; - hHistsxy[name] = filterDF.Histo2D({beamspotName,beamspotName,400,-maxPipeRadius,maxPipeRadius,400,-maxPipeRadius,maxPipeRadius},"xposf","yposf"); - hHistsxyZoom[name] = filterDF.Histo2D({xyname,xyname,100,xMinZoom,xMaxZoom,100,yMinZoom,yMaxZoom},"xposf","yposf"); - hHistsxpx[name] = filterDF.Histo2D({xname,xname,400,xmin,xmax,400,pxmin,pxmax},"xposf","xmomf"); - hHistsypy[name] = filterDF.Histo2D({yname,yname,400,ymin,ymax,400,pymin,pymax},"yposf","ymomf"); - - //Parameters of the pipe - pipeRadii[name] = filterDF.Max("pipeRadiusf").GetValue(); - pipeXPos[name] = filterDF.Max("xdetf").GetValue(); - pipeZPos[name] = filterDF.Max("zdetf").GetValue(); - pipeRotation[name] = filterDF.Max("rotationf").GetValue(); - - //Fit gaussian to the x, y, px and py distributions - auto xhist = hHistsxy[name]->ProjectionX(); - auto yhist = hHistsxy[name]->ProjectionY(); - auto pxhist = hHistsxpx[name]->ProjectionY(); - auto pyhist = hHistsypy[name]->ProjectionY(); - xhist->Fit("gaus","Q"); - yhist->Fit("gaus","Q"); - pxhist->Fit("gaus","Q"); - pyhist->Fit("gaus","Q"); - //Get the fit parameters and errors - xMeanFit[name] = xhist->GetFunction("gaus")->GetParameter(1); - yMeanFit[name] = yhist->GetFunction("gaus")->GetParameter(1); - xMeanFitErr[name] = xhist->GetFunction("gaus")->GetParError(1); - yMeanFitErr[name] = yhist->GetFunction("gaus")->GetParError(1); - xStdDevFit[name] = xhist->GetFunction("gaus")->GetParameter(2); - yStdDevFit[name] = yhist->GetFunction("gaus")->GetParameter(2); - xStdDevFitErr[name] = xhist->GetFunction("gaus")->GetParError(2); - yStdDevFitErr[name] = yhist->GetFunction("gaus")->GetParError(2); - pxMeanFit[name] = pxhist->GetFunction("gaus")->GetParameter(1); - pyMeanFit[name] = pyhist->GetFunction("gaus")->GetParameter(1); - pxMeanFitErr[name] = pxhist->GetFunction("gaus")->GetParError(1); - pyMeanFitErr[name] = pyhist->GetFunction("gaus")->GetParError(1); - pxStdDevFit[name] = pxhist->GetFunction("gaus")->GetParameter(2); - pyStdDevFit[name] = pyhist->GetFunction("gaus")->GetParameter(2); - pxStdDevFitErr[name] = pxhist->GetFunction("gaus")->GetParError(2); - pyStdDevFitErr[name] = pyhist->GetFunction("gaus")->GetParError(2); - - } - - // Create histograms of the beamspot - TCanvas *cXY = new TCanvas("beamspot_canvas","beamspot_canvas",3000,1600); - cXY->Divide(4,2); - int i=1; - for(auto [name,h] : hHistsxy){ - // Get the pipe radius for this histogram - auto pipeRadius = pipeRadii[name]; - 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"); - - // Add zoomed version in the top-right corner - TPad *pad = new TPad("zoomPad", "Zoomed View", 0.65, 0.65, 1.0, 1.0); - pad->SetFillStyle(0); // Transparent background - pad->Draw(); - pad->cd(); - // Draw the zoomed histogram without its title or axis titles - hHistsxyZoom[name]->SetTitle(""); - hHistsxyZoom[name]->GetXaxis()->SetTitle(""); - hHistsxyZoom[name]->GetYaxis()->SetTitle(""); - hHistsxyZoom[name]->Draw("col"); - cXY->cd(); // Return to the main canvas - } - - // x-px canvas - TCanvas *cX = new TCanvas("x_px_canvas","x_px_canvas",3000,1600); - cX->Divide(4,2); - - i=1; - //Write histograms to file - for(auto& h : hHistsxpx){ - cX->cd(i++); - h.second->Draw("col"); - } - - // y-py canvas - TCanvas *cY = new TCanvas("y_py_canvas","y_py_canvas",3000,1600); - cY->Divide(4,2); - - i=1; - for(auto& h : hHistsypy){ - cY->cd(i++); - h.second->Draw("col"); - } - - // Save 2D canvases as pngs - // cXY->SaveAs("beamspot.png"); - // cX->SaveAs("x_px.png"); - // cY->SaveAs("y_py.png"); - - // --------------------------------------------------------------------------- - // Create histograms showing the fitted means and standard deviations of the positions and momenta - // --------------------------------------------------------------------------- - - // Create histograms for fitted X means and standard deviations - TH1F* hFittedXMeans = CreateFittedHistogram("hFittedXMeans", - "Mean X Offset [cm]", - xMeanFit, - xMeanFitErr, - "Pipe ID"); - - TH1F* hFittedXStdDevs = CreateFittedHistogram("hFittedXStdDevs", - "Std Deviation X Offset [cm]", - xStdDevFit, - xStdDevFitErr, - "Pipe ID"); - - // Create histograms for fitted Y means and standard deviations - TH1F* hFittedYMeans = CreateFittedHistogram("hFittedYMeans", - "Mean Y Offset [cm]", - yMeanFit, - yMeanFitErr, - "Pipe ID"); - - TH1F* hFittedYStdDevs = CreateFittedHistogram("hFittedYStdDevs", - "Std Deviation Y Offset [cm]", - yStdDevFit, - yStdDevFitErr, - "Pipe ID"); - - TH1F* hFittedPxMeans = CreateFittedHistogram("hFittedPxMeans", - "Mean Px", - pxMeanFit, - pxMeanFitErr, - "Pipe ID"); - TH1F* hFittedPyMeans = CreateFittedHistogram("hFittedPyMeans", - "Mean Py", - pyMeanFit, - pyMeanFitErr, - "Pipe ID"); - TH1F* hFittedPxStdDevs = CreateFittedHistogram("hFittedPxStdDevs", - "Std Deviation Px", - pxStdDevFit, - pxStdDevFitErr, - "Pipe ID"); - TH1F* hFittedPyStdDevs = CreateFittedHistogram("hFittedPyStdDevs", - "Std Deviation Py", - pyStdDevFit, - pyStdDevFitErr, - "Pipe ID"); - - - // Create a canvas for the fitted beamspot means and standard deviations - TCanvas *cFittedMeans = new TCanvas("cFittedMeans", "Fitted Beamspot Means and Std Deviation", 1200, 800); - cFittedMeans->Divide(2, 2); - cFittedMeans->cd(1); - hFittedXMeans->Draw("E1"); // "E1" draws error bars - cFittedMeans->cd(2); - hFittedXStdDevs->Draw("E1"); // "E1" draws error bars - cFittedMeans->cd(3); - hFittedYMeans->Draw("E1"); // "E1" draws error bars - cFittedMeans->cd(4); - hFittedYStdDevs->Draw("E1"); // "E1" draws error bars - cFittedMeans->SetGrid(); - cFittedMeans->Update(); - // Save the canvas as a PNG file - // cFittedMeans->SaveAs("fitted_means_stddevs.png"); - - // Create a canvas for the fitted momentum means and standard deviations - TCanvas *cFittedMomentumMeans = new TCanvas("cFittedMomentumMeans", "Fitted Momentum Means and Std Deviation", 1200, 800); - cFittedMomentumMeans->Divide(2, 2); - cFittedMomentumMeans->cd(1); - hFittedPxMeans->Draw("E1"); // "E1" draws error bars - cFittedMomentumMeans->cd(2); - hFittedPxStdDevs->Draw("E1"); // "E1" draws error bars - cFittedMomentumMeans->cd(3); - hFittedPyMeans->Draw("E1"); // "E1" draws error bars - cFittedMomentumMeans->cd(4); - hFittedPyStdDevs->Draw("E1"); // "E1" draws error bars - cFittedMomentumMeans->SetGrid(); - cFittedMomentumMeans->Update(); - // Save the canvas as a PNG file - // cFittedMomentumMeans->SaveAs("fitted_momentum_means_stddevs.png"); - - - // ----------------------------------------------------------------------------- - // Create histograms of the beampipe parameters - // ----------------------------------------------------------------------------- - TH1F* hPipeRadii = CreateFittedHistogram("hPipeRadii", - "Radius [cm]", - pipeRadii, - {}, - "Pipe ID"); - TH1F* hPipeXPos = CreateFittedHistogram("hPipeXPos", - "X Position [cm]", - pipeXPos, - {}, - "Pipe ID"); - TH1F* hPipeZPos = CreateFittedHistogram("hPipeZPos", - "Z Position [cm]", - pipeZPos, - {}, - "Pipe ID"); - TH1F* hPipeRotations = CreateFittedHistogram("hPipeRotations", - "Rotation [rad]", - pipeRotation, - {}, - "Pipe ID"); - - // Create a canvas for the pipe parameters - TCanvas *cPipeParams = new TCanvas("cPipeParams", "Pipe Parameters", 1200, 400); - cPipeParams->Divide(4, 1); - cPipeParams->cd(1); - hPipeRadii->Draw(""); - cPipeParams->cd(2); - hPipeXPos->Draw(""); - cPipeParams->cd(3); - hPipeZPos->Draw(""); - cPipeParams->cd(4); - hPipeRotations->Draw(""); - cPipeParams->SetGrid(); - cPipeParams->Update(); - // Save the canvas as a PNG file - // cPipeParams->SaveAs("pipe_parameters.png"); - - - TFile *f = new TFile(outFile,"RECREATE"); - cXY->Write(); - cX->Write(); - cY->Write(); - cFittedMeans->Write(); - cFittedMomentumMeans->Write(); - cPipeParams->Write(); - - f->Close(); - - std::cout << "Saving events to file" << std::endl; - - ROOT::RDF::RSnapshotOptions opts; - // opts.fMode = "UPDATE"; - // d1.Snapshot("events",outFile,{"pipeID","endID","pipeRadius","xpos","ypos","zpos","xmom","ymom","zmom","xpos_global","ypos_global","zpos_global","px_global","py_global","pz_global"},opts); -} diff --git a/benchmarks/beamline/phasespaceGPS.mac b/benchmarks/beamline/phasespaceGPS.mac deleted file mode 100644 index 5d33ceff..00000000 --- a/benchmarks/beamline/phasespaceGPS.mac +++ /dev/null @@ -1,11 +0,0 @@ -/gps/ang/type beam2d -/gps/ang/sigma_x 0.011 rad # 11 mrad -/gps/ang/sigma_y 0.011 rad # 11 mrad -/gps/pos/type Beam -/gps/pos/sigma_x 0.119 mm -/gps/pos/sigma_y 0.0107 mm -/gps/energy 17.846 GeV -#/gps/energy 18 GeV -/gps/particle e- - -/run/beamOn 1000000 \ No newline at end of file From 99633c27222be46530e0f79e4a32c25a18034561 Mon Sep 17 00:00:00 2001 From: simonge Date: Wed, 4 Jun 2025 18:38:36 +0100 Subject: [PATCH 09/36] Prepare for CI testing --- benchmarks/beamline/Snakefile | 51 ++++++++++++++++++++------ benchmarks/beamline/beamlineAnalysis.C | 21 +++++++---- 2 files changed, 53 insertions(+), 19 deletions(-) diff --git a/benchmarks/beamline/Snakefile b/benchmarks/beamline/Snakefile index 88de43c1..6c5639f3 100644 --- a/benchmarks/beamline/Snakefile +++ b/benchmarks/beamline/Snakefile @@ -1,10 +1,8 @@ rule beamline_sim: input: - macro="{test}GPS.mac", + macro="benchmarks/beamline/beamlineGPS.mac", output: - "/scratch/EIC/G4out/beamline/{test}Test{tag}.edm4hep.root", - log: - "/scratch/EIC/G4out/beamline/{test}Test{tag}.edm4hep.root.log", + "sim_output/beamline/beamlineTest{CAMPAIGN}.edm4hep.root", shell: """ exec npsim \ @@ -16,17 +14,46 @@ rule beamline_sim: """ rule beamline_analysis: - input: - data="/scratch/EIC/G4out/beamline/{test}Test{tag}.edm4hep.root", - script="{test}Analysis.C", - output: - "/scratch/EIC/ReconOut/beamline/{test}TestAnalysis{tag}.root" - log: - "/scratch/EIC/ReconOut/beamline/{test}TestAnalysis{tag}.root.log" params: xml=os.getenv("DETECTOR_PATH")+"/epic_ip6_extended.xml", + input: + script="benchmarks/beamline/beamlineAnalysis.C", + data="sim_output/beamline/beamlineTest{CAMPAIGN}.edm4hep.root", + output: + rootfile="sim_output/beamline/analysis/beamlineTestAnalysis{CAMPAIGN}.root", + beamspot_canvas="sim_output/beamline/analysis/beamspot_{CAMPAIGN}.png", + x_px_canvas="sim_output/beamline/analysis/x_px_{CAMPAIGN}.png", + y_py_canvas="sim_output/beamline/analysis/y_py_{CAMPAIGN}.png", + fitted_position_means_stdevs_canvas="sim_output/beamline/analysis/fitted_position_means_stdevs_{CAMPAIGN}.png", + fitted_momentum_means_stdevs_canvas="sim_output/beamline/analysis/fitted_momentum_means_stdevs_{CAMPAIGN}.png", + pipe_parameter_canvas="sim_output/beamline/analysis/pipe_parameter_{CAMPAIGN}.png", shell: """ - root -l -b -q 'analysis.C("{input.data}", "{output}", "{params.xml}")' + root -l -b -q '{input.script}("{input.data}", "{output.rootfile}", "{params.xml}", + "{output.beamspot_canvas}", "{output.x_px_canvas}", "{output.y_py_canvas}", + "{output.fitted_position_means_stdevs_canvas}", "{output.fitted_momentum_means_stdevs_canvas}", + "{output.pipe_parameter_canvas}")' """ +rule beamline: + input: + [ + "sim_output/beamline/analysis/beamlineTestAnalysis{CAMPAIGN}.root", + "sim_output/beamline/analysis/beamspot_{CAMPAIGN}.png", + "sim_output/beamline/analysis/x_px_{CAMPAIGN}.png", + "sim_output/beamline/analysis/y_py_{CAMPAIGN}.png", + "sim_output/beamline/analysis/fitted_position_means_stdevs_{CAMPAIGN}.png", + "sim_output/beamline/analysis/fitted_momentum_means_stdevs_{CAMPAIGN}.png", + "sim_output/beamline/analysis/pipe_parameter_{CAMPAIGN}.png", + ], + output: + directory("results/beamline/{CAMPAIGN}/") + shell: + """ + mkdir {output} + cp {input} {output} + """ + +rule beamline_local: + input: + "results/beamline/local/" \ No newline at end of file diff --git a/benchmarks/beamline/beamlineAnalysis.C b/benchmarks/beamline/beamlineAnalysis.C index 32cd51cd..64b251b2 100644 --- a/benchmarks/beamline/beamlineAnalysis.C +++ b/benchmarks/beamline/beamlineAnalysis.C @@ -21,7 +21,14 @@ using RNode = ROOT::RDF::RNode; void beamlineAnalysis( TString inFile = "/scratch/EIC/G4out/beamline/beamlineTest.edm4hep.root", TString outFile = "output.root", - std::string compactName = "/home/simong/EIC/epic/install/share/epic/epic_ip6_extended.xml"){ + std::string compactName = "/home/simong/EIC/epic/install/share/epic/epic_ip6_extended.xml", + TString beamspotCanvasName = "beamspot_canvas.png", + TString x_pxCanvasName = "x_px_canvas.png", + TString y_pyCanvasName = "y_py_canvas.png", + TString fittedPositionMeansCanvasName = "fitted_means_stddevs.png", + TString fittedMomentumMeansCanvasName = "fitted_momentum_means_stddevs.png", + TString pipeParamsCanvasName = "pipe_parameters.png" + ){ //Set ROOT style gStyle->SetPadLeftMargin(0.1); // Set left margin @@ -305,9 +312,9 @@ void beamlineAnalysis( TString inFile = "/scratch/EIC/G4out/beamline/b } // Save 2D canvases as pngs - cXY->SaveAs("beamspot.png"); - cX->SaveAs("x_px.png"); - cY->SaveAs("y_py.png"); + cXY->SaveAs(beamspotCanvasName); + cX->SaveAs(x_pxCanvasName); + cY->SaveAs(y_pyCanvasName); // --------------------------------------------------------------------------- // Create histograms showing the fitted means and standard deviations of the positions and momenta @@ -375,7 +382,7 @@ void beamlineAnalysis( TString inFile = "/scratch/EIC/G4out/beamline/b cFittedMeans->SetGrid(); cFittedMeans->Update(); // Save the canvas as a PNG file - cFittedMeans->SaveAs("fitted_means_stddevs.png"); + cFittedMeans->SaveAs(fittedPositionMeansCanvasName); // Create a canvas for the fitted momentum means and standard deviations TCanvas *cFittedMomentumMeans = new TCanvas("cFittedMomentumMeans", "Fitted Momentum Means and Std Deviation", 1200, 800); @@ -391,7 +398,7 @@ void beamlineAnalysis( TString inFile = "/scratch/EIC/G4out/beamline/b cFittedMomentumMeans->SetGrid(); cFittedMomentumMeans->Update(); // Save the canvas as a PNG file - cFittedMomentumMeans->SaveAs("fitted_momentum_means_stddevs.png"); + cFittedMomentumMeans->SaveAs(fittedMomentumMeansCanvasName); // ----------------------------------------------------------------------------- @@ -432,7 +439,7 @@ void beamlineAnalysis( TString inFile = "/scratch/EIC/G4out/beamline/b cPipeParams->SetGrid(); cPipeParams->Update(); // Save the canvas as a PNG file - cPipeParams->SaveAs("pipe_parameters.png"); + cPipeParams->SaveAs(pipeParamsCanvasName); TFile *f = new TFile(outFile,"RECREATE"); From 0deb7016325b30c36c7e7f4953d172943525f5e1 Mon Sep 17 00:00:00 2001 From: simonge Date: Wed, 4 Jun 2025 18:38:50 +0100 Subject: [PATCH 10/36] Remove other benchmarks --- .gitlab-ci.yml | 106 +++++++++++++++++++++++++------------------------ Snakefile | 84 +++++++++++++++++++-------------------- 2 files changed, 96 insertions(+), 94 deletions(-) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index c0b49fe0..2b6f3656 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -124,62 +124,64 @@ get_data: - runner_system_failure include: - - local: 'benchmarks/backgrounds/config.yml' - - local: 'benchmarks/backwards_ecal/config.yml' - - local: 'benchmarks/calo_pid/config.yml' - - local: 'benchmarks/ecal_gaps/config.yml' - - local: 'benchmarks/tracking_detectors/config.yml' - - local: 'benchmarks/tracking_performances/config.yml' - - local: 'benchmarks/tracking_performances_dis/config.yml' - - local: 'benchmarks/barrel_ecal/config.yml' - - local: 'benchmarks/barrel_hcal/config.yml' - - local: 'benchmarks/lfhcal/config.yml' - - local: 'benchmarks/zdc/config.yml' - - local: 'benchmarks/zdc_lyso/config.yml' - - local: 'benchmarks/zdc_neutron/config.yml' - - local: 'benchmarks/zdc_photon/config.yml' - - local: 'benchmarks/zdc_pi0/config.yml' - - local: 'benchmarks/material_scan/config.yml' - - local: 'benchmarks/pid/config.yml' - - local: 'benchmarks/timing/config.yml' - - local: 'benchmarks/b0_tracker/config.yml' - - local: 'benchmarks/insert_muon/config.yml' - - local: 'benchmarks/insert_tau/config.yml' - - local: 'benchmarks/zdc_sigma/config.yml' - - local: 'benchmarks/zdc_lambda/config.yml' - - local: 'benchmarks/insert_neutron/config.yml' - - local: 'benchmarks/femc_electron/config.yml' - - local: 'benchmarks/femc_photon/config.yml' - - local: 'benchmarks/femc_pi0/config.yml' + # - local: 'benchmarks/backgrounds/config.yml' + # - local: 'benchmarks/backwards_ecal/config.yml' + # - local: 'benchmarks/calo_pid/config.yml' + # - local: 'benchmarks/ecal_gaps/config.yml' + # - local: 'benchmarks/tracking_detectors/config.yml' + # - local: 'benchmarks/tracking_performances/config.yml' + # - local: 'benchmarks/tracking_performances_dis/config.yml' + # - local: 'benchmarks/barrel_ecal/config.yml' + # - local: 'benchmarks/barrel_hcal/config.yml' + # - local: 'benchmarks/lfhcal/config.yml' + # - local: 'benchmarks/zdc/config.yml' + # - local: 'benchmarks/zdc_lyso/config.yml' + # - local: 'benchmarks/zdc_neutron/config.yml' + # - local: 'benchmarks/zdc_photon/config.yml' + # - local: 'benchmarks/zdc_pi0/config.yml' + # - local: 'benchmarks/material_scan/config.yml' + # - local: 'benchmarks/pid/config.yml' + # - local: 'benchmarks/timing/config.yml' + # - local: 'benchmarks/b0_tracker/config.yml' + # - local: 'benchmarks/insert_muon/config.yml' + # - local: 'benchmarks/insert_tau/config.yml' + # - local: 'benchmarks/zdc_sigma/config.yml' + # - local: 'benchmarks/zdc_lambda/config.yml' + # - local: 'benchmarks/insert_neutron/config.yml' + # - local: 'benchmarks/femc_electron/config.yml' + # - local: 'benchmarks/femc_photon/config.yml' + # - local: 'benchmarks/femc_pi0/config.yml' + - local: 'benchmarks/beamline/config.yml' deploy_results: allow_failure: true stage: deploy needs: - - "collect_results:backgrounds" - - "collect_results:backwards_ecal" - - "collect_results:barrel_ecal" - - "collect_results:barrel_hcal" - - "collect_results:calo_pid" - - "collect_results:ecal_gaps" - - "collect_results:lfhcal" - - "collect_results:material_scan" - - "collect_results:pid" - - "collect_results:tracking_performance" - - "collect_results:tracking_performance_campaigns" - - "collect_results:zdc_sigma" - - "collect_results:zdc_lambda" - - "collect_results:insert_neutron" - - "collect_results:tracking_performances_dis" - - "collect_results:zdc" - - "collect_results:zdc_lyso" - - "collect_results:zdc_neutron" - - "collect_results:insert_muon" - - "collect_results:insert_tau" - - "collect_results:zdc_photon" - - "collect_results:zdc_pi0" - - "collect_results:femc_electron" - - "collect_results:femc_photon" - - "collect_results:femc_pi0" + # - "collect_results:backgrounds" + # - "collect_results:backwards_ecal" + # - "collect_results:barrel_ecal" + # - "collect_results:barrel_hcal" + # - "collect_results:calo_pid" + # - "collect_results:ecal_gaps" + # - "collect_results:lfhcal" + # - "collect_results:material_scan" + # - "collect_results:pid" + # - "collect_results:tracking_performance" + # - "collect_results:tracking_performance_campaigns" + # - "collect_results:zdc_sigma" + # - "collect_results:zdc_lambda" + # - "collect_results:insert_neutron" + # - "collect_results:tracking_performances_dis" + # - "collect_results:zdc" + # - "collect_results:zdc_lyso" + # - "collect_results:zdc_neutron" + # - "collect_results:insert_muon" + # - "collect_results:insert_tau" + # - "collect_results:zdc_photon" + # - "collect_results:zdc_pi0" + # - "collect_results:femc_electron" + # - "collect_results:femc_photon" + # - "collect_results:femc_pi0" + - "collect_results/beamline" script: - snakemake $SNAKEMAKE_FLAGS --cores 1 results/metadata.json - find results -print | sort | tee summary.txt diff --git a/Snakefile b/Snakefile index c001f8b0..3350fb35 100644 --- a/Snakefile +++ b/Snakefile @@ -30,27 +30,27 @@ def find_epic_libraries(): return libs -include: "benchmarks/backgrounds/Snakefile" -include: "benchmarks/backwards_ecal/Snakefile" -include: "benchmarks/barrel_ecal/Snakefile" -include: "benchmarks/calo_pid/Snakefile" -include: "benchmarks/ecal_gaps/Snakefile" -include: "benchmarks/material_scan/Snakefile" -include: "benchmarks/tracking_performances/Snakefile" -include: "benchmarks/tracking_performances_dis/Snakefile" -include: "benchmarks/lfhcal/Snakefile" -include: "benchmarks/zdc_lyso/Snakefile" -include: "benchmarks/zdc_neutron/Snakefile" -include: "benchmarks/insert_muon/Snakefile" -include: "benchmarks/zdc_lambda/Snakefile" -include: "benchmarks/zdc_photon/Snakefile" -include: "benchmarks/zdc_pi0/Snakefile" -include: "benchmarks/zdc_sigma/Snakefile" -include: "benchmarks/insert_neutron/Snakefile" -include: "benchmarks/insert_tau/Snakefile" -include: "benchmarks/femc_electron/Snakefile" -include: "benchmarks/femc_photon/Snakefile" -include: "benchmarks/femc_pi0/Snakefile" +#include: "benchmarks/backgrounds/Snakefile" +#include: "benchmarks/backwards_ecal/Snakefile" +#include: "benchmarks/barrel_ecal/Snakefile" +#include: "benchmarks/calo_pid/Snakefile" +#include: "benchmarks/ecal_gaps/Snakefile" +#include: "benchmarks/material_scan/Snakefile" +#include: "benchmarks/tracking_performances/Snakefile" +#include: "benchmarks/tracking_performances_dis/Snakefile" +#include: "benchmarks/lfhcal/Snakefile" +#include: "benchmarks/zdc_lyso/Snakefile" +#include: "benchmarks/zdc_neutron/Snakefile" +#include: "benchmarks/insert_muon/Snakefile" +#include: "benchmarks/zdc_lambda/Snakefile" +#include: "benchmarks/zdc_photon/Snakefile" +#include: "benchmarks/zdc_pi0/Snakefile" +#include: "benchmarks/zdc_sigma/Snakefile" +#include: "benchmarks/insert_neutron/Snakefile" +#include: "benchmarks/insert_tau/Snakefile" +#include: "benchmarks/femc_electron/Snakefile" +#include: "benchmarks/femc_photon/Snakefile" +#include: "benchmarks/femc_pi0/Snakefile" include: "benchmarks/beamline/Snakefile" use_s3 = config["remote_provider"].lower() == "s3" @@ -121,24 +121,24 @@ awk -f {input.converter} {input.notebook} > {output} """ -rule metadata: - output: - "results/metadata.json" - shell: - """ -cat > {output} < {output} < Date: Wed, 4 Jun 2025 19:33:55 +0100 Subject: [PATCH 11/36] Add missing config.yml --- benchmarks/beamline/config.yml | 26 ++++++++++++++++++++++++++ 1 file changed, 26 insertions(+) create mode 100644 benchmarks/beamline/config.yml diff --git a/benchmarks/beamline/config.yml b/benchmarks/beamline/config.yml new file mode 100644 index 00000000..996a5fb4 --- /dev/null +++ b/benchmarks/beamline/config.yml @@ -0,0 +1,26 @@ +sim:beamline: + extends: .det_benchmark + stage: simulate + script: + - | + snakemake --cache --cores 5 \ + sim_output/beamline/beamlineTestlocal.edm4hep.root + +bench:beamline: + extends: .det_benchmark + stage: benchmarks + needs: + - ["sim:beamline"] + script: + - snakemake $SNAKEMAKE_FLAGS --cores 3 beamline_local + +collect_results:beamline: + extends: .det_benchmark + stage: collect + needs: + - "bench:beamline" + script: + - ls -lrht + - mv results{,_save}/ # move results directory out of the way to preserve it + - snakemake $SNAKEMAKE_FLAGS --cores 1 --delete-all-output beamline_local + - mv results{_save,}/ From 7a583553acbc7aa5f819d26b7b8aa1cac9edd1ed Mon Sep 17 00:00:00 2001 From: simonge Date: Wed, 4 Jun 2025 19:39:19 +0100 Subject: [PATCH 12/36] Correct typo --- .gitlab-ci.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 2b6f3656..bd6d0383 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -181,7 +181,7 @@ deploy_results: # - "collect_results:femc_electron" # - "collect_results:femc_photon" # - "collect_results:femc_pi0" - - "collect_results/beamline" + - "collect_results:beamline" script: - snakemake $SNAKEMAKE_FLAGS --cores 1 results/metadata.json - find results -print | sort | tee summary.txt From 5cd62be1c807b637816ce19a98a939dcaa3f0a0c Mon Sep 17 00:00:00 2001 From: simonge Date: Tue, 17 Jun 2025 17:16:26 +0100 Subject: [PATCH 13/36] Re-enable other benchmarks and update cores --- Snakefile | 42 +++++++++++++++++----------------- benchmarks/beamline/config.yml | 4 ++-- 2 files changed, 23 insertions(+), 23 deletions(-) diff --git a/Snakefile b/Snakefile index 3350fb35..5e374afd 100644 --- a/Snakefile +++ b/Snakefile @@ -30,27 +30,27 @@ def find_epic_libraries(): return libs -#include: "benchmarks/backgrounds/Snakefile" -#include: "benchmarks/backwards_ecal/Snakefile" -#include: "benchmarks/barrel_ecal/Snakefile" -#include: "benchmarks/calo_pid/Snakefile" -#include: "benchmarks/ecal_gaps/Snakefile" -#include: "benchmarks/material_scan/Snakefile" -#include: "benchmarks/tracking_performances/Snakefile" -#include: "benchmarks/tracking_performances_dis/Snakefile" -#include: "benchmarks/lfhcal/Snakefile" -#include: "benchmarks/zdc_lyso/Snakefile" -#include: "benchmarks/zdc_neutron/Snakefile" -#include: "benchmarks/insert_muon/Snakefile" -#include: "benchmarks/zdc_lambda/Snakefile" -#include: "benchmarks/zdc_photon/Snakefile" -#include: "benchmarks/zdc_pi0/Snakefile" -#include: "benchmarks/zdc_sigma/Snakefile" -#include: "benchmarks/insert_neutron/Snakefile" -#include: "benchmarks/insert_tau/Snakefile" -#include: "benchmarks/femc_electron/Snakefile" -#include: "benchmarks/femc_photon/Snakefile" -#include: "benchmarks/femc_pi0/Snakefile" +include: "benchmarks/backgrounds/Snakefile" +include: "benchmarks/backwards_ecal/Snakefile" +include: "benchmarks/barrel_ecal/Snakefile" +include: "benchmarks/calo_pid/Snakefile" +include: "benchmarks/ecal_gaps/Snakefile" +include: "benchmarks/material_scan/Snakefile" +include: "benchmarks/tracking_performances/Snakefile" +include: "benchmarks/tracking_performances_dis/Snakefile" +include: "benchmarks/lfhcal/Snakefile" +include: "benchmarks/zdc_lyso/Snakefile" +include: "benchmarks/zdc_neutron/Snakefile" +include: "benchmarks/insert_muon/Snakefile" +include: "benchmarks/zdc_lambda/Snakefile" +include: "benchmarks/zdc_photon/Snakefile" +include: "benchmarks/zdc_pi0/Snakefile" +include: "benchmarks/zdc_sigma/Snakefile" +include: "benchmarks/insert_neutron/Snakefile" +include: "benchmarks/insert_tau/Snakefile" +include: "benchmarks/femc_electron/Snakefile" +include: "benchmarks/femc_photon/Snakefile" +include: "benchmarks/femc_pi0/Snakefile" include: "benchmarks/beamline/Snakefile" use_s3 = config["remote_provider"].lower() == "s3" diff --git a/benchmarks/beamline/config.yml b/benchmarks/beamline/config.yml index 996a5fb4..ce147524 100644 --- a/benchmarks/beamline/config.yml +++ b/benchmarks/beamline/config.yml @@ -3,7 +3,7 @@ sim:beamline: stage: simulate script: - | - snakemake --cache --cores 5 \ + snakemake --cache --cores 2 \ sim_output/beamline/beamlineTestlocal.edm4hep.root bench:beamline: @@ -12,7 +12,7 @@ bench:beamline: needs: - ["sim:beamline"] script: - - snakemake $SNAKEMAKE_FLAGS --cores 3 beamline_local + - snakemake $SNAKEMAKE_FLAGS --cores 2 beamline_local collect_results:beamline: extends: .det_benchmark From 312ca425c4b70607cb4091137f4d7bb40956a55c Mon Sep 17 00:00:00 2001 From: simonge Date: Wed, 18 Jun 2025 18:29:42 +0100 Subject: [PATCH 14/36] Add some pass fail checks --- benchmarks/beamline/beamlineAnalysis.C | 51 ++++++++++++++++++++++++-- benchmarks/beamline/beamlineGPS.mac | 6 +-- 2 files changed, 50 insertions(+), 7 deletions(-) diff --git a/benchmarks/beamline/beamlineAnalysis.C b/benchmarks/beamline/beamlineAnalysis.C index 64b251b2..851c4fd0 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; -void beamlineAnalysis( TString inFile = "/scratch/EIC/G4out/beamline/beamlineTest.edm4hep.root", +int beamlineAnalysis( TString inFile = "/scratch/EIC/G4out/beamline/beamlineTest.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", @@ -46,6 +46,8 @@ void beamlineAnalysis( TString inFile = "/scratch/EIC/G4out/beamline/b //Set implicit multi-threading ROOT::EnableImplicitMT(); + int pass = 0; + //Load the detector config dd4hep::Detector& detector = dd4hep::Detector::getInstance(); detector.fromCompact(compactName); @@ -56,6 +58,9 @@ void beamlineAnalysis( TString inFile = "/scratch/EIC/G4out/beamline/b //Set number of entries to process // d1 = d1.Range(0,1000); + int nEntries = d1.Count().GetValue(); + float acceptableLoss = 0.99; // Set the acceptable loss percentage to 1% + float acceptableEntries = nEntries * acceptableLoss; //Get the collection std::string readoutName = "BackwardsBeamlineHits"; @@ -120,7 +125,7 @@ void beamlineAnalysis( TString inFile = "/scratch/EIC/G4out/beamline/b } else{ std::cout << "Collection " << readoutName << " not found in file" << std::endl; - return; + return 1; } // Calculate the maximum pipe radius for plotting @@ -199,6 +204,32 @@ void beamlineAnalysis( TString inFile = "/scratch/EIC/G4out/beamline/b auto pxmaxf = filterDF.Max("xmomf").GetValue(); auto pyminf = filterDF.Min("ymomf").GetValue(); auto pymaxf = filterDF.Max("ymomf").GetValue(); + // If the min and max values are equal, set them to a small value + if(xminf == xmaxf) { + xminf -= 0.01; + xmaxf += 0.01; + pass = 1; // Set pass to 1 if xminf == xmaxf + std::cout << "Warning: xminf == xmaxf for pipe ID " << i << ", likely no hits." << std::endl; + } + if(yminf == ymaxf) { + yminf -= 0.01; + ymaxf += 0.01; + pass = 1; // Set pass to 1 if yminf == ymaxf + std::cout << "Warning: yminf == ymaxf for pipe ID " << i << ", likely no hits." << std::endl; + } + if(pxminf == pxmaxf) { + pxminf -= 0.01; + pxmaxf += 0.01; + pass = 1; // Set pass to 1 if pxminf == pxmaxf + std::cout << "Warning: pxminf == pxmaxf for pipe ID " << i << ", likely no hits." << std::endl; + } + if(pyminf == pymaxf) { + pyminf -= 0.01; + pymaxf += 0.01; + pass = 1; // Set pass to 1 if pyminf == pymaxf + std::cout << "Warning: pyminf == pymaxf for pipe ID " << i << ", likely no hits." << std::endl; + } + // Calculate means and standard deviations xMeans[name] = filterDF.Mean("xposf").GetValue(); yMeans[name] = filterDF.Mean("yposf").GetValue(); @@ -227,6 +258,7 @@ void beamlineAnalysis( TString inFile = "/scratch/EIC/G4out/beamline/b //Parameters of the pipe pipeRadii[name] = filterDF.Max("pipeRadiusf").GetValue(); + std::cout << "Pipe ID: " << name << " Radius: " << pipeRadii[name] << std::endl; pipeXPos[name] = filterDF.Max("xdetf").GetValue(); pipeZPos[name] = filterDF.Max("zdetf").GetValue(); pipeRotation[name] = filterDF.Max("rotationf").GetValue(); @@ -265,6 +297,14 @@ void beamlineAnalysis( TString inFile = "/scratch/EIC/G4out/beamline/b cXY->Divide(4,2); int i=1; for(auto [name,h] : hHistsxy){ + + // Check histogram entries + if(h->GetEntries() < acceptableEntries){ + std::cout << "Warning: Only " << h->GetEntries()/nEntries << " of particles contributing to histogram " << name + << " , which is below the accepted threshold of " << acceptableEntries << std::endl; + pass = 1; + } + // Get the pipe radius for this histogram auto pipeRadius = pipeRadii[name]; cXY->cd(i++); @@ -452,9 +492,12 @@ void beamlineAnalysis( TString inFile = "/scratch/EIC/G4out/beamline/b f->Close(); - std::cout << "Saving events to file" << std::endl; + std::cout << "Analysis complete. Results saved to " << outFile << std::endl; + return pass; + + // std::cout << "Saving events to file" << std::endl; - ROOT::RDF::RSnapshotOptions opts; + // ROOT::RDF::RSnapshotOptions opts; // opts.fMode = "UPDATE"; // d1.Snapshot("events",outFile,{"pipeID","endID","pipeRadius","xpos","ypos","zpos","xmom","ymom","zmom","xpos_global","ypos_global","zpos_global","px_global","py_global","pz_global"},opts); } diff --git a/benchmarks/beamline/beamlineGPS.mac b/benchmarks/beamline/beamlineGPS.mac index 8ecd6bbb..d3cdc21b 100644 --- a/benchmarks/beamline/beamlineGPS.mac +++ b/benchmarks/beamline/beamlineGPS.mac @@ -4,8 +4,8 @@ /gps/pos/type Beam /gps/pos/sigma_x 0.119 mm /gps/pos/sigma_y 0.0107 mm -/gps/energy 17.846 GeV -#/gps/energy 18 GeV +#/gps/energy 17.846 GeV +/gps/energy 18 GeV /gps/particle e- -/run/beamOn 100000 \ No newline at end of file +/run/beamOn 10000 \ No newline at end of file From 02e0c0e03af4972c8c0327de57c0ebf9268ff6c2 Mon Sep 17 00:00:00 2001 From: simonge Date: Wed, 18 Jun 2025 21:14:07 +0100 Subject: [PATCH 15/36] Set sim and analysis dir by variable --- benchmarks/beamline/Snakefile | 39 +++++++++++++++++++---------------- 1 file changed, 21 insertions(+), 18 deletions(-) diff --git a/benchmarks/beamline/Snakefile b/benchmarks/beamline/Snakefile index 6c5639f3..5722a943 100644 --- a/benchmarks/beamline/Snakefile +++ b/benchmarks/beamline/Snakefile @@ -1,8 +1,11 @@ -rule beamline_sim: +SIMOUTDIR="sim_output/beamline/" +ANALYSISDIR=SIMOUTDIR+"analysis/" + +rule beamline_steering_sim: input: macro="benchmarks/beamline/beamlineGPS.mac", output: - "sim_output/beamline/beamlineTest{CAMPAIGN}.edm4hep.root", + SIMOUTDIR+"beamlineTest{CAMPAIGN}.edm4hep.root", shell: """ exec npsim \ @@ -13,20 +16,20 @@ rule beamline_sim: --outputFile {output} \ """ -rule beamline_analysis: +rule beamline_steering_analysis: params: xml=os.getenv("DETECTOR_PATH")+"/epic_ip6_extended.xml", input: script="benchmarks/beamline/beamlineAnalysis.C", - data="sim_output/beamline/beamlineTest{CAMPAIGN}.edm4hep.root", + data=SIMOUTDIR+"beamlineTest{CAMPAIGN}.edm4hep.root", output: - rootfile="sim_output/beamline/analysis/beamlineTestAnalysis{CAMPAIGN}.root", - beamspot_canvas="sim_output/beamline/analysis/beamspot_{CAMPAIGN}.png", - x_px_canvas="sim_output/beamline/analysis/x_px_{CAMPAIGN}.png", - y_py_canvas="sim_output/beamline/analysis/y_py_{CAMPAIGN}.png", - fitted_position_means_stdevs_canvas="sim_output/beamline/analysis/fitted_position_means_stdevs_{CAMPAIGN}.png", - fitted_momentum_means_stdevs_canvas="sim_output/beamline/analysis/fitted_momentum_means_stdevs_{CAMPAIGN}.png", - pipe_parameter_canvas="sim_output/beamline/analysis/pipe_parameter_{CAMPAIGN}.png", + rootfile=ANALYSISDIR+"beamlineTestAnalysis{CAMPAIGN}.root", + beamspot_canvas=ANALYSISDIR+"beamspot_{CAMPAIGN}.png", + x_px_canvas=ANALYSISDIR+"x_px_{CAMPAIGN}.png", + y_py_canvas=ANALYSISDIR+"y_py_{CAMPAIGN}.png", + fitted_position_means_stdevs_canvas=ANALYSISDIR+"fitted_position_means_stdevs_{CAMPAIGN}.png", + fitted_momentum_means_stdevs_canvas=ANALYSISDIR+"fitted_momentum_means_stdevs_{CAMPAIGN}.png", + pipe_parameter_canvas=ANALYSISDIR+"pipe_parameter_{CAMPAIGN}.png", shell: """ root -l -b -q '{input.script}("{input.data}", "{output.rootfile}", "{params.xml}", @@ -38,13 +41,13 @@ rule beamline_analysis: rule beamline: input: [ - "sim_output/beamline/analysis/beamlineTestAnalysis{CAMPAIGN}.root", - "sim_output/beamline/analysis/beamspot_{CAMPAIGN}.png", - "sim_output/beamline/analysis/x_px_{CAMPAIGN}.png", - "sim_output/beamline/analysis/y_py_{CAMPAIGN}.png", - "sim_output/beamline/analysis/fitted_position_means_stdevs_{CAMPAIGN}.png", - "sim_output/beamline/analysis/fitted_momentum_means_stdevs_{CAMPAIGN}.png", - "sim_output/beamline/analysis/pipe_parameter_{CAMPAIGN}.png", + ANALYSISDIR+"beamlineTestAnalysis{CAMPAIGN}.root", + ANALYSISDIR+"beamspot_{CAMPAIGN}.png", + ANALYSISDIR+"x_px_{CAMPAIGN}.png", + ANALYSISDIR+"y_py_{CAMPAIGN}.png", + ANALYSISDIR+"fitted_position_means_stdevs_{CAMPAIGN}.png", + ANALYSISDIR+"fitted_momentum_means_stdevs_{CAMPAIGN}.png", + ANALYSISDIR+"pipe_parameter_{CAMPAIGN}.png", ], output: directory("results/beamline/{CAMPAIGN}/") From e39d0f9d4ef4f1271a9b904bce5c75428614b33c Mon Sep 17 00:00:00 2001 From: simonge Date: Thu, 19 Jun 2025 10:38:28 +0100 Subject: [PATCH 16/36] Revert exclusion of other benchmarks --- .gitlab-ci.yml | 104 ++++++++++++++++++++++++------------------------- Snakefile | 42 ++++++++++---------- 2 files changed, 73 insertions(+), 73 deletions(-) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index bd6d0383..31456e7e 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -124,64 +124,64 @@ get_data: - runner_system_failure include: - # - local: 'benchmarks/backgrounds/config.yml' - # - local: 'benchmarks/backwards_ecal/config.yml' - # - local: 'benchmarks/calo_pid/config.yml' - # - local: 'benchmarks/ecal_gaps/config.yml' - # - local: 'benchmarks/tracking_detectors/config.yml' - # - local: 'benchmarks/tracking_performances/config.yml' - # - local: 'benchmarks/tracking_performances_dis/config.yml' - # - local: 'benchmarks/barrel_ecal/config.yml' - # - local: 'benchmarks/barrel_hcal/config.yml' - # - local: 'benchmarks/lfhcal/config.yml' - # - local: 'benchmarks/zdc/config.yml' - # - local: 'benchmarks/zdc_lyso/config.yml' - # - local: 'benchmarks/zdc_neutron/config.yml' - # - local: 'benchmarks/zdc_photon/config.yml' - # - local: 'benchmarks/zdc_pi0/config.yml' - # - local: 'benchmarks/material_scan/config.yml' - # - local: 'benchmarks/pid/config.yml' - # - local: 'benchmarks/timing/config.yml' - # - local: 'benchmarks/b0_tracker/config.yml' - # - local: 'benchmarks/insert_muon/config.yml' - # - local: 'benchmarks/insert_tau/config.yml' - # - local: 'benchmarks/zdc_sigma/config.yml' - # - local: 'benchmarks/zdc_lambda/config.yml' - # - local: 'benchmarks/insert_neutron/config.yml' - # - local: 'benchmarks/femc_electron/config.yml' - # - local: 'benchmarks/femc_photon/config.yml' - # - local: 'benchmarks/femc_pi0/config.yml' + - local: 'benchmarks/backgrounds/config.yml' + - local: 'benchmarks/backwards_ecal/config.yml' - local: 'benchmarks/beamline/config.yml' + - local: 'benchmarks/calo_pid/config.yml' + - local: 'benchmarks/ecal_gaps/config.yml' + - local: 'benchmarks/tracking_detectors/config.yml' + - local: 'benchmarks/tracking_performances/config.yml' + - local: 'benchmarks/tracking_performances_dis/config.yml' + - local: 'benchmarks/barrel_ecal/config.yml' + - local: 'benchmarks/barrel_hcal/config.yml' + - local: 'benchmarks/lfhcal/config.yml' + - local: 'benchmarks/zdc/config.yml' + - local: 'benchmarks/zdc_lyso/config.yml' + - local: 'benchmarks/zdc_neutron/config.yml' + - local: 'benchmarks/zdc_photon/config.yml' + - local: 'benchmarks/zdc_pi0/config.yml' + - local: 'benchmarks/material_scan/config.yml' + - local: 'benchmarks/pid/config.yml' + - local: 'benchmarks/timing/config.yml' + - local: 'benchmarks/b0_tracker/config.yml' + - local: 'benchmarks/insert_muon/config.yml' + - local: 'benchmarks/insert_tau/config.yml' + - local: 'benchmarks/zdc_sigma/config.yml' + - local: 'benchmarks/zdc_lambda/config.yml' + - local: 'benchmarks/insert_neutron/config.yml' + - local: 'benchmarks/femc_electron/config.yml' + - local: 'benchmarks/femc_photon/config.yml' + - local: 'benchmarks/femc_pi0/config.yml' deploy_results: allow_failure: true stage: deploy needs: - # - "collect_results:backgrounds" - # - "collect_results:backwards_ecal" - # - "collect_results:barrel_ecal" - # - "collect_results:barrel_hcal" - # - "collect_results:calo_pid" - # - "collect_results:ecal_gaps" - # - "collect_results:lfhcal" - # - "collect_results:material_scan" - # - "collect_results:pid" - # - "collect_results:tracking_performance" - # - "collect_results:tracking_performance_campaigns" - # - "collect_results:zdc_sigma" - # - "collect_results:zdc_lambda" - # - "collect_results:insert_neutron" - # - "collect_results:tracking_performances_dis" - # - "collect_results:zdc" - # - "collect_results:zdc_lyso" - # - "collect_results:zdc_neutron" - # - "collect_results:insert_muon" - # - "collect_results:insert_tau" - # - "collect_results:zdc_photon" - # - "collect_results:zdc_pi0" - # - "collect_results:femc_electron" - # - "collect_results:femc_photon" - # - "collect_results:femc_pi0" + - "collect_results:backgrounds" + - "collect_results:backwards_ecal" + - "collect_results:barrel_ecal" + - "collect_results:barrel_hcal" - "collect_results:beamline" + - "collect_results:calo_pid" + - "collect_results:ecal_gaps" + - "collect_results:lfhcal" + - "collect_results:material_scan" + - "collect_results:pid" + - "collect_results:tracking_performance" + - "collect_results:tracking_performance_campaigns" + - "collect_results:zdc_sigma" + - "collect_results:zdc_lambda" + - "collect_results:insert_neutron" + - "collect_results:tracking_performances_dis" + - "collect_results:zdc" + - "collect_results:zdc_lyso" + - "collect_results:zdc_neutron" + - "collect_results:insert_muon" + - "collect_results:insert_tau" + - "collect_results:zdc_photon" + - "collect_results:zdc_pi0" + - "collect_results:femc_electron" + - "collect_results:femc_photon" + - "collect_results:femc_pi0" script: - snakemake $SNAKEMAKE_FLAGS --cores 1 results/metadata.json - find results -print | sort | tee summary.txt diff --git a/Snakefile b/Snakefile index 5e374afd..d70cd9e4 100644 --- a/Snakefile +++ b/Snakefile @@ -121,24 +121,24 @@ awk -f {input.converter} {input.notebook} > {output} """ -#rule metadata: -# output: -# "results/metadata.json" -# shell: -# """ -#cat > {output} < {output} < Date: Thu, 19 Jun 2025 11:33:55 +0100 Subject: [PATCH 17/36] Review suggestions --- Snakefile | 2 +- benchmarks/beamline/Snakefile | 16 +++++++--------- benchmarks/beamline/config.yml | 2 +- 3 files changed, 9 insertions(+), 11 deletions(-) diff --git a/Snakefile b/Snakefile index d70cd9e4..4e856e1d 100644 --- a/Snakefile +++ b/Snakefile @@ -33,6 +33,7 @@ def find_epic_libraries(): include: "benchmarks/backgrounds/Snakefile" include: "benchmarks/backwards_ecal/Snakefile" include: "benchmarks/barrel_ecal/Snakefile" +include: "benchmarks/beamline/Snakefile" include: "benchmarks/calo_pid/Snakefile" include: "benchmarks/ecal_gaps/Snakefile" include: "benchmarks/material_scan/Snakefile" @@ -51,7 +52,6 @@ include: "benchmarks/insert_tau/Snakefile" include: "benchmarks/femc_electron/Snakefile" include: "benchmarks/femc_photon/Snakefile" include: "benchmarks/femc_pi0/Snakefile" -include: "benchmarks/beamline/Snakefile" use_s3 = config["remote_provider"].lower() == "s3" use_xrootd = config["remote_provider"].lower() == "xrootd" diff --git a/benchmarks/beamline/Snakefile b/benchmarks/beamline/Snakefile index 5722a943..bfba8748 100644 --- a/benchmarks/beamline/Snakefile +++ b/benchmarks/beamline/Snakefile @@ -40,15 +40,13 @@ rule beamline_steering_analysis: rule beamline: input: - [ - ANALYSISDIR+"beamlineTestAnalysis{CAMPAIGN}.root", - ANALYSISDIR+"beamspot_{CAMPAIGN}.png", - ANALYSISDIR+"x_px_{CAMPAIGN}.png", - ANALYSISDIR+"y_py_{CAMPAIGN}.png", - ANALYSISDIR+"fitted_position_means_stdevs_{CAMPAIGN}.png", - ANALYSISDIR+"fitted_momentum_means_stdevs_{CAMPAIGN}.png", - ANALYSISDIR+"pipe_parameter_{CAMPAIGN}.png", - ], + ANALYSISDIR+"beamlineTestAnalysis{CAMPAIGN}.root", + ANALYSISDIR+"beamspot_{CAMPAIGN}.png", + ANALYSISDIR+"x_px_{CAMPAIGN}.png", + ANALYSISDIR+"y_py_{CAMPAIGN}.png", + ANALYSISDIR+"fitted_position_means_stdevs_{CAMPAIGN}.png", + ANALYSISDIR+"fitted_momentum_means_stdevs_{CAMPAIGN}.png", + ANALYSISDIR+"pipe_parameter_{CAMPAIGN}.png", output: directory("results/beamline/{CAMPAIGN}/") shell: diff --git a/benchmarks/beamline/config.yml b/benchmarks/beamline/config.yml index ce147524..ef6b84e4 100644 --- a/benchmarks/beamline/config.yml +++ b/benchmarks/beamline/config.yml @@ -3,7 +3,7 @@ sim:beamline: stage: simulate script: - | - snakemake --cache --cores 2 \ + snakemake $SNAKEMAKE_FLAGS --cores 2 \ sim_output/beamline/beamlineTestlocal.edm4hep.root bench:beamline: From 81fad99121202b5349dd84191e5ba06d3e2454cc Mon Sep 17 00:00:00 2001 From: Dmitry Kalinkin Date: Thu, 19 Jun 2025 18:45:56 -0400 Subject: [PATCH 18/36] Snakefile: fix for out-of-tree running --- benchmarks/beamline/Snakefile | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/benchmarks/beamline/Snakefile b/benchmarks/beamline/Snakefile index bfba8748..a43b5054 100644 --- a/benchmarks/beamline/Snakefile +++ b/benchmarks/beamline/Snakefile @@ -3,7 +3,7 @@ ANALYSISDIR=SIMOUTDIR+"analysis/" rule beamline_steering_sim: input: - macro="benchmarks/beamline/beamlineGPS.mac", + macro=workflow.source_path("beamlineGPS.mac"), output: SIMOUTDIR+"beamlineTest{CAMPAIGN}.edm4hep.root", shell: @@ -20,7 +20,7 @@ rule beamline_steering_analysis: params: xml=os.getenv("DETECTOR_PATH")+"/epic_ip6_extended.xml", input: - script="benchmarks/beamline/beamlineAnalysis.C", + script=workflow.source_path("beamlineAnalysis.C"), data=SIMOUTDIR+"beamlineTest{CAMPAIGN}.edm4hep.root", output: rootfile=ANALYSISDIR+"beamlineTestAnalysis{CAMPAIGN}.root", @@ -57,4 +57,4 @@ rule beamline: rule beamline_local: input: - "results/beamline/local/" \ No newline at end of file + "results/beamline/local/" From 6b21ff11333ae22cf9c4dc380b5213d34653abe0 Mon Sep 17 00:00:00 2001 From: Dmitry Kalinkin Date: Thu, 19 Jun 2025 18:46:43 -0400 Subject: [PATCH 19/36] Snakefile restore indentation --- Snakefile | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/Snakefile b/Snakefile index 4e856e1d..2277a081 100644 --- a/Snakefile +++ b/Snakefile @@ -122,21 +122,21 @@ awk -f {input.converter} {input.notebook} > {output} rule metadata: - output: - "results/metadata.json" - shell: - """ + output: + "results/metadata.json" + shell: + """ cat > {output} < Date: Mon, 23 Jun 2025 17:19:59 +0100 Subject: [PATCH 20/36] Add header to inputs too --- benchmarks/beamline/Snakefile | 1 + 1 file changed, 1 insertion(+) diff --git a/benchmarks/beamline/Snakefile b/benchmarks/beamline/Snakefile index a43b5054..b9478bc6 100644 --- a/benchmarks/beamline/Snakefile +++ b/benchmarks/beamline/Snakefile @@ -21,6 +21,7 @@ rule beamline_steering_analysis: xml=os.getenv("DETECTOR_PATH")+"/epic_ip6_extended.xml", input: script=workflow.source_path("beamlineAnalysis.C"), + header=workflow.source_path("shared_functions.h"), data=SIMOUTDIR+"beamlineTest{CAMPAIGN}.edm4hep.root", output: rootfile=ANALYSISDIR+"beamlineTestAnalysis{CAMPAIGN}.root", From c7e268f10cf7c59fe51cb3a5aaa3634a4474606d Mon Sep 17 00:00:00 2001 From: simonge Date: Tue, 20 May 2025 13:51:19 +0100 Subject: [PATCH 21/36] Add low-q2 phase space electron tests --- benchmarks/beamline/phasespaceAnalysis.C | 253 +++++++++++++++++++++++ benchmarks/beamline/phasespaceGPS.mac | 13 ++ 2 files changed, 266 insertions(+) create mode 100644 benchmarks/beamline/phasespaceAnalysis.C create mode 100644 benchmarks/beamline/phasespaceGPS.mac diff --git a/benchmarks/beamline/phasespaceAnalysis.C b/benchmarks/beamline/phasespaceAnalysis.C new file mode 100644 index 00000000..4792e653 --- /dev/null +++ b/benchmarks/beamline/phasespaceAnalysis.C @@ -0,0 +1,253 @@ +// Script to plot the x and y positions and momentum of beamline particles as they pass through the magnets + +#include +#include "DD4hep/Detector.h" +#include "DDRec/CellIDPositionConverter.h" +#include "edm4hep/MCParticleCollection.h" +#include "edm4hep/SimTrackerHitCollection.h" +#include "edm4hep/SimCalorimeterHitCollection.h" +#include "edm4hep/Vector3f.h" +#include "edm4hep/Vector3d.h" +#include "ROOT/RDataFrame.hxx" +#include "ROOT/RDF/RInterface.hxx" +#include "ROOT/RVec.hxx" +#include "shared_functions.h" +#include "TCanvas.h" +#include "TStyle.h" + + +using RVecS = ROOT::VecOps::RVec; +using RNode = ROOT::RDF::RNode; + +void phasespaceAnalysis( TString inFile = "/scratch/EIC/G4out/beamline/beamlineTest.edm4hep.root", + TString outFile = "output.root", + std::string compactName = "/home/simong/EIC/epic/install/share/epic/epic_ip6_extended.xml"){ + + //Set ROOT style + gStyle->SetPadLeftMargin(0.1); // Set left margin + gStyle->SetPadRightMargin(0.0); // Set right margin + gStyle->SetPadTopMargin(0.0); // Set top margin + gStyle->SetPadBottomMargin(0.1);// Set bottom margin + gStyle->SetTitleAlign(13); + gStyle->SetTitleX(0.12); // Place the title on the top right + gStyle->SetTitleY(0.985); // Place the title on the top right + gStyle->SetTitleSize(0.08, "t"); + gStyle->SetTitleXOffset(1.0); // Adjust y-axis title offset + gStyle->SetTitleYOffset(1.0); // Adjust y-axis title offset + gStyle->SetOptStat(0); + + //Set implicit multi-threading + ROOT::EnableImplicitMT(); + + //Load the detector config + dd4hep::Detector& detector = dd4hep::Detector::getInstance(); + detector.fromCompact(compactName); + + ROOT::RDataFrame d0("events",inFile, {"BackwardsBeamlineHits"}); + RNode d1 = d0; + RVecS colNames = d1.GetColumnNames(); + + //Set number of entries to process + // d1 = d1.Range(0,1000); + + //Get the collection + std::string readoutName = "BackwardsBeamlineHits"; + + std::cout << "Running lazy RDataframe execution" << std::endl; + + if(Any(colNames==readoutName)){ + + d1 = d1.Define("pipeID",getSubID("pipe",detector),{readoutName}) + .Define("endID",getSubID("end",detector),{readoutName}) + .Define("pipeParameters",getVolumeParametersFromCellID(detector),{readoutName}) + .Define("pipeRadius",[](const ROOT::VecOps::RVec& params) { + ROOT::VecOps::RVec radii; + for (const auto& param : params) { + radii.push_back(param.radius); + } + return radii; + }, {"pipeParameters"}) + .Define("xdet",[](const ROOT::VecOps::RVec& params) { + ROOT::VecOps::RVec xPos; + for (const auto& param : params) { + xPos.push_back(param.xPos); + } + return xPos; + }, {"pipeParameters"}) + .Define("zdet",[](const ROOT::VecOps::RVec& params) { + ROOT::VecOps::RVec zPos; + for (const auto& param : params) { + zPos.push_back(param.zPos); + } + return zPos; + }, {"pipeParameters"}) + .Define("rotation",[](const ROOT::VecOps::RVec& params) { + ROOT::VecOps::RVec rotation; + for (const auto& param : params) { + rotation.push_back(param.rotation); + } + return rotation; + }, {"pipeParameters"}); + + + //global x,y,z position and momentum + d1 = d1.Define("xpos_global","BackwardsBeamlineHits.position.x") + .Define("ypos_global","BackwardsBeamlineHits.position.y") + .Define("zpos_global","BackwardsBeamlineHits.position.z") + .Define("px_global","BackwardsBeamlineHits.momentum.x") + .Define("py_global","BackwardsBeamlineHits.momentum.y") + .Define("pz_global","BackwardsBeamlineHits.momentum.z"); + + d1 = d1.Define("hitPosMom",globalToLocal(detector),{readoutName}) + .Define("xpos","hitPosMom[0]") + .Define("ypos","hitPosMom[1]") + .Define("zpos","hitPosMom[2]") + .Define("xmomMag","hitPosMom[3]") + .Define("ymomMag","hitPosMom[4]") + .Define("zmomMag","hitPosMom[5]") + .Define("momMag","sqrt(xmomMag*xmomMag+ymomMag*ymomMag+zmomMag*zmomMag)") + .Define("xmom","xmomMag/momMag") + .Define("ymom","ymomMag/momMag") + .Define("zmom","zmomMag/momMag"); + + } + else{ + std::cout << "Collection " << readoutName << " not found in file" << std::endl; + return; + } + + // Calculate the maximum pipe radius for plotting + auto maxPipeRadius = 1.2*d1.Max("pipeRadius").GetValue(); + + std::cout << "Executing Analysis and creating histograms" << std::endl; + + //Create array of histogram results + std::map> hHistsxy; + std::map> hHistsxyZoom; + std::map> hHistsxpx; + std::map> hHistsypy; + + std::map xMeans; + std::map yMeans; + std::map xStdDevs; + std::map yStdDevs; + std::map pxMeans; + std::map pyMeans; + std::map pxStdDevs; + std::map pyStdDevs; + + //Fit paremeter and error maps + std::map xMeanFit; + std::map yMeanFit; + std::map xMeanFitErr; + std::map yMeanFitErr; + std::map xStdDevFit; + std::map yStdDevFit; + std::map xStdDevFitErr; + std::map yStdDevFitErr; + std::map pxMeanFit; + std::map pyMeanFit; + std::map pxMeanFitErr; + std::map pyMeanFitErr; + std::map pxStdDevFit; + std::map pyStdDevFit; + std::map pxStdDevFitErr; + std::map pyStdDevFitErr; + + std::map pipeRadii; + std::map pipeXPos; + std::map pipeZPos; + std::map pipeRotation; + + auto xmin = d1.Min("xpos").GetValue(); + auto xmax = d1.Max("xpos").GetValue(); + auto ymin = d1.Min("ypos").GetValue(); + auto ymax = d1.Max("ypos").GetValue(); + auto pxmin = d1.Min("xmom").GetValue(); + auto pxmax = d1.Max("xmom").GetValue(); + auto pymin = d1.Min("ymom").GetValue(); + auto pymax = d1.Max("ymom").GetValue(); + + //Create histograms + for(int i=0; i<=7; i++){ + + std::string name = "pipeID"; + name += std::to_string(i); + auto filterDF = d1.Define("xposf","xpos["+std::to_string(i)+"]") + .Define("yposf","ypos["+std::to_string(i)+"]") + .Define("xmomf","xmom["+std::to_string(i)+"]") + .Define("ymomf","ymom["+std::to_string(i)+"]") + .Define("pipeRadiusf","pipeRadius["+std::to_string(i)+"]") + .Define("xdetf","xdet["+std::to_string(i)+"]") + .Define("zdetf","zdet["+std::to_string(i)+"]") + .Define("rotationf","rotation["+std::to_string(i)+"]"); + + + //Calculate Min and Max values + auto xminf = filterDF.Min("xposf").GetValue(); + auto xmaxf = filterDF.Max("xposf").GetValue(); + auto yminf = filterDF.Min("yposf").GetValue(); + auto ymaxf = filterDF.Max("yposf").GetValue(); + auto pxminf = filterDF.Min("xmomf").GetValue(); + auto pxmaxf = filterDF.Max("xmomf").GetValue(); + auto pyminf = filterDF.Min("ymomf").GetValue(); + auto pymaxf = filterDF.Max("ymomf").GetValue(); + // Calculate means and standard deviations + xMeans[name] = filterDF.Mean("xposf").GetValue(); + yMeans[name] = filterDF.Mean("yposf").GetValue(); + xStdDevs[name] = filterDF.StdDev("xposf").GetValue(); + yStdDevs[name] = filterDF.StdDev("yposf").GetValue(); + pxMeans[name] = filterDF.Mean("xmomf").GetValue(); + pyMeans[name] = filterDF.Mean("ymomf").GetValue(); + pxStdDevs[name] = filterDF.StdDev("xmomf").GetValue(); + pyStdDevs[name] = filterDF.StdDev("ymomf").GetValue(); + + + TString beamspotName = "Beamspot ID"+std::to_string(i)+";x offset [cm]; y offset [cm]"; + TString xyname = name+";x Offset [cm]; y Offset [cm]"; + TString xname = name+";x Offset [cm]; x trajectory component"; + TString yname = name+";y Offset [cm]; y trajectory component"; + hHistsxy[name] = filterDF.Histo2D({beamspotName,beamspotName,400,-maxPipeRadius,maxPipeRadius,400,-maxPipeRadius,maxPipeRadius},"xposf","yposf"); + + //Parameters of the pipe + pipeRadii[name] = filterDF.Max("pipeRadiusf").GetValue(); + pipeXPos[name] = filterDF.Max("xdetf").GetValue(); + pipeZPos[name] = filterDF.Max("zdetf").GetValue(); + pipeRotation[name] = filterDF.Max("rotationf").GetValue(); + + } + + // Create histograms of the beamspot + TCanvas *cXY = new TCanvas("beamspot_canvas","beamspot_canvas",3000,1600); + cXY->Divide(4,2); + int i=1; + for(auto [name,h] : hHistsxy){ + // Get the pipe radius for this histogram + auto pipeRadius = pipeRadii[name]; + 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"); + + } + + // Save 2D canvases as pngs + cXY->SaveAs("phasespace_in_beampipe.png"); + + + TFile *f = new TFile(outFile,"RECREATE"); + cXY->Write(); + + f->Close(); + + std::cout << "Saving events to file" << std::endl; + + ROOT::RDF::RSnapshotOptions opts; + opts.fMode = "UPDATE"; + d1.Snapshot("events",outFile,{"pipeID","endID","pipeRadius","xpos","ypos","zpos","xmom","ymom","zmom","xpos_global","ypos_global","zpos_global","px_global","py_global","pz_global"},opts); +} diff --git a/benchmarks/beamline/phasespaceGPS.mac b/benchmarks/beamline/phasespaceGPS.mac new file mode 100644 index 00000000..81b1a814 --- /dev/null +++ b/benchmarks/beamline/phasespaceGPS.mac @@ -0,0 +1,13 @@ +/gps/ang/type beam2d +/gps/ang/sigma_x 0.0011 rad # 11 mrad +/gps/ang/sigma_y 0.0011 rad # 11 mrad +/gps/pos/type Beam +/gps/pos/sigma_x 0.119 mm +/gps/pos/sigma_y 0.0107 mm +#/gps/energy 17.846 GeV +/gps/ene/type Lin +/gps/ene/min 6 GeV +/gps/ene/max 17 GeV +/gps/particle e- + +/run/beamOn 100000 \ No newline at end of file From c964d40ea6c2225c695222ff1b9db4fe5e0d0f90 Mon Sep 17 00:00:00 2001 From: simonge Date: Thu, 19 Jun 2025 10:36:37 +0100 Subject: [PATCH 22/36] Add acceptance sim --- benchmarks/beamline/Snakefile | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/benchmarks/beamline/Snakefile b/benchmarks/beamline/Snakefile index b9478bc6..750a8884 100644 --- a/benchmarks/beamline/Snakefile +++ b/benchmarks/beamline/Snakefile @@ -16,6 +16,21 @@ rule beamline_steering_sim: --outputFile {output} \ """ +rule beamline_acceptance_sim: + input: + macro="benchmarks/beamline/phasespaceGPS.mac", + output: + SIMOUTDIR+"phasespaceTest{CAMPAIGN}.edm4hep.root", + shell: + """ + exec npsim \ + --runType run \ + --enableG4GPS \ + --macroFile {input.macro} \ + --compactFile $DETECTOR_PATH/epic_ip6_extended.xml \ + --outputFile {output} \ + """ + rule beamline_steering_analysis: params: xml=os.getenv("DETECTOR_PATH")+"/epic_ip6_extended.xml", From 8b98ed6f700ae37fce2f5d0cb7ce19be4656c139 Mon Sep 17 00:00:00 2001 From: simonge Date: Fri, 20 Jun 2025 14:07:51 +0100 Subject: [PATCH 23/36] Make simulation run with correct range and 1000x faster --- benchmarks/beamline/Snakefile | 24 +++++++++++++++++++++++- benchmarks/beamline/phasespaceGPS.mac | 13 +++++++------ 2 files changed, 30 insertions(+), 7 deletions(-) diff --git a/benchmarks/beamline/Snakefile b/benchmarks/beamline/Snakefile index 750a8884..ba1b8507 100644 --- a/benchmarks/beamline/Snakefile +++ b/benchmarks/beamline/Snakefile @@ -14,6 +14,7 @@ rule beamline_steering_sim: --macroFile {input.macro} \ --compactFile $DETECTOR_PATH/epic_ip6_extended.xml \ --outputFile {output} \ + --physics.rangecut 1000000 \ """ rule beamline_acceptance_sim: @@ -29,6 +30,7 @@ rule beamline_acceptance_sim: --macroFile {input.macro} \ --compactFile $DETECTOR_PATH/epic_ip6_extended.xml \ --outputFile {output} \ + --physics.rangecut 1000000 \ """ rule beamline_steering_analysis: @@ -54,6 +56,22 @@ rule beamline_steering_analysis: "{output.pipe_parameter_canvas}")' """ +rule beamline_acceptance_analysis: + params: + xml=os.getenv("DETECTOR_PATH")+"/epic_ip6_extended.xml", + input: + script="benchmarks/beamline/phasespaceAnalysis.C", + data=SIMOUTDIR+"phasespaceTest{CAMPAIGN}.edm4hep.root", + output: + rootfile=ANALYSISDIR+"phasespaceTestAnalysis{CAMPAIGN}.root", + beampipe_canvas=ANALYSISDIR+"phasespace_in_beampipe_{CAMPAIGN}.png", + etheta_canvas=ANALYSISDIR+"phasespace_energy_theta_{CAMPAIGN}.png", + etheta_acceptance_canvas=ANALYSISDIR+"phasespace_energy_theta_acceptance_{CAMPAIGN}.png", + shell: + """ + root -l -b -q '{input.script}("{input.data}", "{output.rootfile}", "{params.xml}", "{output.beampipe_canvas}","{output.etheta_canvas}","{output.etheta_acceptance_canvas}")' + """ + rule beamline: input: ANALYSISDIR+"beamlineTestAnalysis{CAMPAIGN}.root", @@ -62,7 +80,11 @@ rule beamline: ANALYSISDIR+"y_py_{CAMPAIGN}.png", ANALYSISDIR+"fitted_position_means_stdevs_{CAMPAIGN}.png", ANALYSISDIR+"fitted_momentum_means_stdevs_{CAMPAIGN}.png", - ANALYSISDIR+"pipe_parameter_{CAMPAIGN}.png", + ANALYSISDIR+"pipe_parameter_{CAMPAIGN}.png", + ANALYSISDIR+"phasespaceTestAnalysis{CAMPAIGN}.root", + ANALYSISDIR+"phasespace_in_beampipe_{CAMPAIGN}.png", + ANALYSISDIR+"phasespace_energy_theta_{CAMPAIGN}.png", + ANALYSISDIR+"phasespace_energy_theta_acceptance_{CAMPAIGN}.png", output: directory("results/beamline/{CAMPAIGN}/") shell: diff --git a/benchmarks/beamline/phasespaceGPS.mac b/benchmarks/beamline/phasespaceGPS.mac index 81b1a814..59b48263 100644 --- a/benchmarks/beamline/phasespaceGPS.mac +++ b/benchmarks/beamline/phasespaceGPS.mac @@ -1,13 +1,14 @@ + /gps/ang/type beam2d -/gps/ang/sigma_x 0.0011 rad # 11 mrad -/gps/ang/sigma_y 0.0011 rad # 11 mrad +/gps/ang/sigma_x 0.011 rad # 11 mrad +/gps/ang/sigma_y 0.011 rad # 11 mrad /gps/pos/type Beam /gps/pos/sigma_x 0.119 mm /gps/pos/sigma_y 0.0107 mm -#/gps/energy 17.846 GeV +/gps/particle e- /gps/ene/type Lin +/gps/ene/gradient 0.00001 /gps/ene/min 6 GeV -/gps/ene/max 17 GeV -/gps/particle e- +/gps/ene/max 18 GeV -/run/beamOn 100000 \ No newline at end of file +/run/beamOn 1000000 \ No newline at end of file From 31fe8ef6efc649492cb55963c2ac18d9dbcc0c47 Mon Sep 17 00:00:00 2001 From: simonge Date: Fri, 20 Jun 2025 14:09:14 +0100 Subject: [PATCH 24/36] Add outputs from script --- benchmarks/beamline/phasespaceAnalysis.C | 139 +++++++++++++++++------ 1 file changed, 106 insertions(+), 33 deletions(-) diff --git a/benchmarks/beamline/phasespaceAnalysis.C b/benchmarks/beamline/phasespaceAnalysis.C index 4792e653..20e74b80 100644 --- a/benchmarks/beamline/phasespaceAnalysis.C +++ b/benchmarks/beamline/phasespaceAnalysis.C @@ -19,9 +19,12 @@ using RVecS = ROOT::VecOps::RVec; using RNode = ROOT::RDF::RNode; -void phasespaceAnalysis( TString inFile = "/scratch/EIC/G4out/beamline/beamlineTest.edm4hep.root", - TString outFile = "output.root", - std::string compactName = "/home/simong/EIC/epic/install/share/epic/epic_ip6_extended.xml"){ +int phasespaceAnalysis( TString inFile = "/scratch/EIC/G4out/beamline/beamlineTest.edm4hep.root", + TString outFile = "output.root", + std::string compactName = "/home/simong/EIC/epic/install/share/epic/epic_ip6_extended.xml", + TString beampipeCanvasName = "phasespace_in_beampipe.png", + TString EThetaCanvasName = "phasespace_energy_theta.png", + TString EThetaAccCanvasName= "phasespace_energy_theta_acceptance.png") { //Set ROOT style gStyle->SetPadLeftMargin(0.1); // Set left margin @@ -36,8 +39,10 @@ void phasespaceAnalysis( TString inFile = "/scratch/EIC/G4out/beamline/bea gStyle->SetTitleYOffset(1.0); // Adjust y-axis title offset gStyle->SetOptStat(0); + int pass = 0; + //Set implicit multi-threading - ROOT::EnableImplicitMT(); + // ROOT::EnableImplicitMT(); //Load the detector config dd4hep::Detector& detector = dd4hep::Detector::getInstance(); @@ -51,10 +56,44 @@ void phasespaceAnalysis( TString inFile = "/scratch/EIC/G4out/beamline/bea // d1 = d1.Range(0,1000); //Get the collection + std::string mcParticlesName = "MCParticles"; std::string readoutName = "BackwardsBeamlineHits"; std::cout << "Running lazy RDataframe execution" << std::endl; + if(Any(colNames==mcParticlesName)){ + + // Get theta and energy of the particles + d1 = d1.Define("theta",[](const vector& mcParticles) { + // Calculate theta from momentum components as angle from negative z axis + double px = mcParticles[0].momentum.x; + double py = mcParticles[0].momentum.y; + double pz = mcParticles[0].momentum.z; + double p = std::sqrt(px*px + py*py + pz*pz); + double theta = M_PI-std::acos(pz / p); // Angle from the z-axis + + return theta; + }, {mcParticlesName}) + .Define("energy",[](const vector& mcParticles) { + + //Calculate energy from mass and momentum + double mass = mcParticles[0].mass; + double px = mcParticles[0].momentum.x; + double py = mcParticles[0].momentum.y; + double pz = mcParticles[0].momentum.z; + double energy = std::sqrt(px*px + py*py + pz*pz + mass*mass); + + return energy; + }, {mcParticlesName}); + + } + else{ + std::cout << "Collection " << mcParticlesName << " not found in file" << std::endl; + return 1; + } + + auto totalETheta = d1.Histo2D({"Energy vs Theta","Energy vs Theta",100,4,18,100,0,0.011},"energy","theta"); + if(Any(colNames==readoutName)){ d1 = d1.Define("pipeID",getSubID("pipe",detector),{readoutName}) @@ -91,13 +130,14 @@ void phasespaceAnalysis( TString inFile = "/scratch/EIC/G4out/beamline/bea //global x,y,z position and momentum - d1 = d1.Define("xpos_global","BackwardsBeamlineHits.position.x") - .Define("ypos_global","BackwardsBeamlineHits.position.y") - .Define("zpos_global","BackwardsBeamlineHits.position.z") - .Define("px_global","BackwardsBeamlineHits.momentum.x") - .Define("py_global","BackwardsBeamlineHits.momentum.y") - .Define("pz_global","BackwardsBeamlineHits.momentum.z"); - + d1 = d1 .Define("NHits","BackwardsBeamlineHits.size()") + .Alias("xpos_global","BackwardsBeamlineHits.position.x") + .Alias("ypos_global","BackwardsBeamlineHits.position.y") + .Alias("zpos_global","BackwardsBeamlineHits.position.z") + .Alias("px_global","BackwardsBeamlineHits.momentum.x") + .Alias("py_global","BackwardsBeamlineHits.momentum.y") + .Alias("pz_global","BackwardsBeamlineHits.momentum.z"); + d1 = d1.Define("hitPosMom",globalToLocal(detector),{readoutName}) .Define("xpos","hitPosMom[0]") .Define("ypos","hitPosMom[1]") @@ -113,7 +153,7 @@ void phasespaceAnalysis( TString inFile = "/scratch/EIC/G4out/beamline/bea } else{ std::cout << "Collection " << readoutName << " not found in file" << std::endl; - return; + return 1; } // Calculate the maximum pipe radius for plotting @@ -123,9 +163,8 @@ void phasespaceAnalysis( TString inFile = "/scratch/EIC/G4out/beamline/bea //Create array of histogram results std::map> hHistsxy; - std::map> hHistsxyZoom; - std::map> hHistsxpx; - std::map> hHistsypy; + std::map> hHistsETheta; + std::map xMeans; std::map yMeans; @@ -172,17 +211,20 @@ void phasespaceAnalysis( TString inFile = "/scratch/EIC/G4out/beamline/bea for(int i=0; i<=7; i++){ std::string name = "pipeID"; - name += std::to_string(i); - auto filterDF = d1.Define("xposf","xpos["+std::to_string(i)+"]") - .Define("yposf","ypos["+std::to_string(i)+"]") - .Define("xmomf","xmom["+std::to_string(i)+"]") - .Define("ymomf","ymom["+std::to_string(i)+"]") - .Define("pipeRadiusf","pipeRadius["+std::to_string(i)+"]") - .Define("xdetf","xdet["+std::to_string(i)+"]") - .Define("zdetf","zdet["+std::to_string(i)+"]") - .Define("rotationf","rotation["+std::to_string(i)+"]"); + std::string str_i = std::to_string(i); + name += str_i; + auto filterDF = d1.Filter(std::to_string(i+1)+"<=NHits" ) + .Define("xposf","xpos["+str_i+"]") + .Define("yposf","ypos["+str_i+"]") + .Define("xmomf","xmom["+str_i+"]") + .Define("ymomf","ymom["+str_i+"]") + .Define("pipeRadiusf","pipeRadius["+str_i+"]") + .Define("xdetf","xdet["+str_i+"]") + .Define("zdetf","zdet["+str_i+"]") + .Define("rotationf","rotation["+str_i+"]"); - + // auto report = filterDF.Display({"NHits","theta","energy"}); + // report->Print(); //Calculate Min and Max values auto xminf = filterDF.Min("xposf").GetValue(); auto xmaxf = filterDF.Max("xposf").GetValue(); @@ -203,14 +245,20 @@ void phasespaceAnalysis( TString inFile = "/scratch/EIC/G4out/beamline/bea pyStdDevs[name] = filterDF.StdDev("ymomf").GetValue(); - TString beamspotName = "Beamspot ID"+std::to_string(i)+";x offset [cm]; y offset [cm]"; + TString beamspotName = "Beamspot ID"+str_i+";x offset [cm]; y offset [cm]"; TString xyname = name+";x Offset [cm]; y Offset [cm]"; TString xname = name+";x Offset [cm]; x trajectory component"; TString yname = name+";y Offset [cm]; y trajectory component"; hHistsxy[name] = filterDF.Histo2D({beamspotName,beamspotName,400,-maxPipeRadius,maxPipeRadius,400,-maxPipeRadius,maxPipeRadius},"xposf","yposf"); + auto extraFilterDF = filterDF.Filter(std::to_string(i+1)+"==NHits" ); + TString EThetaName = "Energy vs Theta ID"+str_i+";Energy [GeV]; Theta [rad]"; + TString nameETheta = name+";Energy [GeV]; Theta [rad]"; + hHistsETheta[name] = extraFilterDF.Histo2D({EThetaName,EThetaName,100,4,18,100,0,0.011},"energy","theta"); + //Parameters of the pipe - pipeRadii[name] = filterDF.Max("pipeRadiusf").GetValue(); + pipeRadii[name] = filterDF.Max("pipeRadiusf").GetValue(); + std::cout << "Pipe ID: " << name << " Radius: " << pipeRadii[name] << " " << filterDF.Min("pipeRadiusf").GetValue() << std::endl; pipeXPos[name] = filterDF.Max("xdetf").GetValue(); pipeZPos[name] = filterDF.Max("zdetf").GetValue(); pipeRotation[name] = filterDF.Max("rotationf").GetValue(); @@ -236,18 +284,43 @@ void phasespaceAnalysis( TString inFile = "/scratch/EIC/G4out/beamline/bea } - // Save 2D canvases as pngs - cXY->SaveAs("phasespace_in_beampipe.png"); + // Cnavas for energy vs theta + TCanvas *cETheta = new TCanvas("energy_theta_canvas","energy_theta_canvas",3000,1600); + cETheta->Divide(4,2); + i=1; + for(auto [name,h] : hHistsETheta){ + cETheta->cd(i++); + h->Draw("colz"); + } + + // Canvas for energy vs theta acceptance + TCanvas *cEThetaAcc = new TCanvas("energy_theta_acceptance_canvas","energy_theta_acceptance_canvas",3000,1600); + cEThetaAcc->Divide(4,2); + i=1; + for(auto [name,h] : hHistsETheta){ + cEThetaAcc->cd(i++); + h->Divide(totalETheta.GetPtr()); + h->Draw("colz"); + } + // Save 2D canvases as pngs + cXY->SaveAs(beampipeCanvasName); + cETheta->SaveAs(EThetaCanvasName); + cEThetaAcc->SaveAs(EThetaAccCanvasName); TFile *f = new TFile(outFile,"RECREATE"); cXY->Write(); + cETheta->Write(); + cEThetaAcc->Write(); f->Close(); - std::cout << "Saving events to file" << std::endl; + std::cout << "Analysis complete. Results saved to " << outFile << std::endl; + return pass; + + // std::cout << "Saving events to file" << std::endl; - ROOT::RDF::RSnapshotOptions opts; - opts.fMode = "UPDATE"; - d1.Snapshot("events",outFile,{"pipeID","endID","pipeRadius","xpos","ypos","zpos","xmom","ymom","zmom","xpos_global","ypos_global","zpos_global","px_global","py_global","pz_global"},opts); + // ROOT::RDF::RSnapshotOptions opts; + // opts.fMode = "UPDATE"; + // d1.Snapshot("events",outFile,{"pipeID","endID","pipeRadius","xpos","ypos","zpos","xmom","ymom","zmom","xpos_global","ypos_global","zpos_global","px_global","py_global","pz_global"},opts); } From ccddf53e19a1b374b9d0a7e3672e7a35ec87313b Mon Sep 17 00:00:00 2001 From: simonge Date: Mon, 23 Jun 2025 22:57:14 +0100 Subject: [PATCH 25/36] Change code inputs to workflow source path --- benchmarks/beamline/Snakefile | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/benchmarks/beamline/Snakefile b/benchmarks/beamline/Snakefile index ba1b8507..90943244 100644 --- a/benchmarks/beamline/Snakefile +++ b/benchmarks/beamline/Snakefile @@ -19,7 +19,7 @@ rule beamline_steering_sim: rule beamline_acceptance_sim: input: - macro="benchmarks/beamline/phasespaceGPS.mac", + macro=workflow.source_path("phasespaceGPS.mac"), output: SIMOUTDIR+"phasespaceTest{CAMPAIGN}.edm4hep.root", shell: @@ -60,7 +60,8 @@ rule beamline_acceptance_analysis: params: xml=os.getenv("DETECTOR_PATH")+"/epic_ip6_extended.xml", input: - script="benchmarks/beamline/phasespaceAnalysis.C", + script=workflow.source_path("phasespaceAnalysis.C"), + header=workflow.source_path("shared_functions.h"), data=SIMOUTDIR+"phasespaceTest{CAMPAIGN}.edm4hep.root", output: rootfile=ANALYSISDIR+"phasespaceTestAnalysis{CAMPAIGN}.root", From f44d6900721abab8064b5aab5fbd54ac78d410b3 Mon Sep 17 00:00:00 2001 From: simonge Date: Mon, 23 Jun 2025 23:02:59 +0100 Subject: [PATCH 26/36] rename phasespace to acceptance --- benchmarks/beamline/Snakefile | 24 +++++++++---------- ...sespaceAnalysis.C => acceptanceAnalysis.C} | 2 +- .../{phasespaceGPS.mac => acceptanceGPS.mac} | 0 3 files changed, 13 insertions(+), 13 deletions(-) rename benchmarks/beamline/{phasespaceAnalysis.C => acceptanceAnalysis.C} (99%) rename benchmarks/beamline/{phasespaceGPS.mac => acceptanceGPS.mac} (100%) diff --git a/benchmarks/beamline/Snakefile b/benchmarks/beamline/Snakefile index 90943244..a1cc0911 100644 --- a/benchmarks/beamline/Snakefile +++ b/benchmarks/beamline/Snakefile @@ -19,9 +19,9 @@ rule beamline_steering_sim: rule beamline_acceptance_sim: input: - macro=workflow.source_path("phasespaceGPS.mac"), + macro=workflow.source_path("acceptanceGPS.mac"), output: - SIMOUTDIR+"phasespaceTest{CAMPAIGN}.edm4hep.root", + SIMOUTDIR+"acceptanceTest{CAMPAIGN}.edm4hep.root", shell: """ exec npsim \ @@ -60,14 +60,14 @@ rule beamline_acceptance_analysis: params: xml=os.getenv("DETECTOR_PATH")+"/epic_ip6_extended.xml", input: - script=workflow.source_path("phasespaceAnalysis.C"), + script=workflow.source_path("acceptanceAnalysis.C"), header=workflow.source_path("shared_functions.h"), - data=SIMOUTDIR+"phasespaceTest{CAMPAIGN}.edm4hep.root", + data=SIMOUTDIR+"acceptanceTest{CAMPAIGN}.edm4hep.root", output: - rootfile=ANALYSISDIR+"phasespaceTestAnalysis{CAMPAIGN}.root", - beampipe_canvas=ANALYSISDIR+"phasespace_in_beampipe_{CAMPAIGN}.png", - etheta_canvas=ANALYSISDIR+"phasespace_energy_theta_{CAMPAIGN}.png", - etheta_acceptance_canvas=ANALYSISDIR+"phasespace_energy_theta_acceptance_{CAMPAIGN}.png", + rootfile=ANALYSISDIR+"acceptanceTestAnalysis{CAMPAIGN}.root", + beampipe_canvas=ANALYSISDIR+"acceptance_in_beampipe_{CAMPAIGN}.png", + etheta_canvas=ANALYSISDIR+"acceptance_energy_theta_{CAMPAIGN}.png", + etheta_acceptance_canvas=ANALYSISDIR+"acceptance_energy_theta_acceptance_{CAMPAIGN}.png", shell: """ root -l -b -q '{input.script}("{input.data}", "{output.rootfile}", "{params.xml}", "{output.beampipe_canvas}","{output.etheta_canvas}","{output.etheta_acceptance_canvas}")' @@ -82,10 +82,10 @@ rule beamline: ANALYSISDIR+"fitted_position_means_stdevs_{CAMPAIGN}.png", ANALYSISDIR+"fitted_momentum_means_stdevs_{CAMPAIGN}.png", ANALYSISDIR+"pipe_parameter_{CAMPAIGN}.png", - ANALYSISDIR+"phasespaceTestAnalysis{CAMPAIGN}.root", - ANALYSISDIR+"phasespace_in_beampipe_{CAMPAIGN}.png", - ANALYSISDIR+"phasespace_energy_theta_{CAMPAIGN}.png", - ANALYSISDIR+"phasespace_energy_theta_acceptance_{CAMPAIGN}.png", + ANALYSISDIR+"acceptanceTestAnalysis{CAMPAIGN}.root", + ANALYSISDIR+"acceptance_in_beampipe_{CAMPAIGN}.png", + ANALYSISDIR+"acceptance_energy_theta_{CAMPAIGN}.png", + ANALYSISDIR+"acceptance_energy_theta_acceptance_{CAMPAIGN}.png", output: directory("results/beamline/{CAMPAIGN}/") shell: diff --git a/benchmarks/beamline/phasespaceAnalysis.C b/benchmarks/beamline/acceptanceAnalysis.C similarity index 99% rename from benchmarks/beamline/phasespaceAnalysis.C rename to benchmarks/beamline/acceptanceAnalysis.C index 20e74b80..7645caec 100644 --- a/benchmarks/beamline/phasespaceAnalysis.C +++ b/benchmarks/beamline/acceptanceAnalysis.C @@ -19,7 +19,7 @@ using RVecS = ROOT::VecOps::RVec; using RNode = ROOT::RDF::RNode; -int phasespaceAnalysis( TString inFile = "/scratch/EIC/G4out/beamline/beamlineTest.edm4hep.root", +int acceptanceAnalysis( TString inFile = "/scratch/EIC/G4out/beamline/acceptanceTest.edm4hep.root", TString outFile = "output.root", std::string compactName = "/home/simong/EIC/epic/install/share/epic/epic_ip6_extended.xml", TString beampipeCanvasName = "phasespace_in_beampipe.png", diff --git a/benchmarks/beamline/phasespaceGPS.mac b/benchmarks/beamline/acceptanceGPS.mac similarity index 100% rename from benchmarks/beamline/phasespaceGPS.mac rename to benchmarks/beamline/acceptanceGPS.mac From 4df2119ba4a9046dd325ea4d48a9e6a414390b77 Mon Sep 17 00:00:00 2001 From: simonge Date: Mon, 23 Jun 2025 23:07:26 +0100 Subject: [PATCH 27/36] Remove unused code --- benchmarks/beamline/acceptanceAnalysis.C | 110 +---------------------- 1 file changed, 4 insertions(+), 106 deletions(-) diff --git a/benchmarks/beamline/acceptanceAnalysis.C b/benchmarks/beamline/acceptanceAnalysis.C index 7645caec..19cb21ff 100644 --- a/benchmarks/beamline/acceptanceAnalysis.C +++ b/benchmarks/beamline/acceptanceAnalysis.C @@ -105,50 +105,16 @@ int acceptanceAnalysis( TString inFile = "/scratch/EIC/G4out/beamlin radii.push_back(param.radius); } return radii; - }, {"pipeParameters"}) - .Define("xdet",[](const ROOT::VecOps::RVec& params) { - ROOT::VecOps::RVec xPos; - for (const auto& param : params) { - xPos.push_back(param.xPos); - } - return xPos; - }, {"pipeParameters"}) - .Define("zdet",[](const ROOT::VecOps::RVec& params) { - ROOT::VecOps::RVec zPos; - for (const auto& param : params) { - zPos.push_back(param.zPos); - } - return zPos; - }, {"pipeParameters"}) - .Define("rotation",[](const ROOT::VecOps::RVec& params) { - ROOT::VecOps::RVec rotation; - for (const auto& param : params) { - rotation.push_back(param.rotation); - } - return rotation; }, {"pipeParameters"}); //global x,y,z position and momentum - d1 = d1 .Define("NHits","BackwardsBeamlineHits.size()") - .Alias("xpos_global","BackwardsBeamlineHits.position.x") - .Alias("ypos_global","BackwardsBeamlineHits.position.y") - .Alias("zpos_global","BackwardsBeamlineHits.position.z") - .Alias("px_global","BackwardsBeamlineHits.momentum.x") - .Alias("py_global","BackwardsBeamlineHits.momentum.y") - .Alias("pz_global","BackwardsBeamlineHits.momentum.z"); + d1 = d1 .Define("NHits","BackwardsBeamlineHits.size()"); d1 = d1.Define("hitPosMom",globalToLocal(detector),{readoutName}) .Define("xpos","hitPosMom[0]") .Define("ypos","hitPosMom[1]") - .Define("zpos","hitPosMom[2]") - .Define("xmomMag","hitPosMom[3]") - .Define("ymomMag","hitPosMom[4]") - .Define("zmomMag","hitPosMom[5]") - .Define("momMag","sqrt(xmomMag*xmomMag+ymomMag*ymomMag+zmomMag*zmomMag)") - .Define("xmom","xmomMag/momMag") - .Define("ymom","ymomMag/momMag") - .Define("zmom","zmomMag/momMag"); + .Define("zpos","hitPosMom[2]"); } else{ @@ -166,46 +132,7 @@ int acceptanceAnalysis( TString inFile = "/scratch/EIC/G4out/beamlin std::map> hHistsETheta; - std::map xMeans; - std::map yMeans; - std::map xStdDevs; - std::map yStdDevs; - std::map pxMeans; - std::map pyMeans; - std::map pxStdDevs; - std::map pyStdDevs; - - //Fit paremeter and error maps - std::map xMeanFit; - std::map yMeanFit; - std::map xMeanFitErr; - std::map yMeanFitErr; - std::map xStdDevFit; - std::map yStdDevFit; - std::map xStdDevFitErr; - std::map yStdDevFitErr; - std::map pxMeanFit; - std::map pyMeanFit; - std::map pxMeanFitErr; - std::map pyMeanFitErr; - std::map pxStdDevFit; - std::map pyStdDevFit; - std::map pxStdDevFitErr; - std::map pyStdDevFitErr; - std::map pipeRadii; - std::map pipeXPos; - std::map pipeZPos; - std::map pipeRotation; - - auto xmin = d1.Min("xpos").GetValue(); - auto xmax = d1.Max("xpos").GetValue(); - auto ymin = d1.Min("ypos").GetValue(); - auto ymax = d1.Max("ypos").GetValue(); - auto pxmin = d1.Min("xmom").GetValue(); - auto pxmax = d1.Max("xmom").GetValue(); - auto pymin = d1.Min("ymom").GetValue(); - auto pymax = d1.Max("ymom").GetValue(); //Create histograms for(int i=0; i<=7; i++){ @@ -216,34 +143,8 @@ int acceptanceAnalysis( TString inFile = "/scratch/EIC/G4out/beamlin auto filterDF = d1.Filter(std::to_string(i+1)+"<=NHits" ) .Define("xposf","xpos["+str_i+"]") .Define("yposf","ypos["+str_i+"]") - .Define("xmomf","xmom["+str_i+"]") - .Define("ymomf","ymom["+str_i+"]") - .Define("pipeRadiusf","pipeRadius["+str_i+"]") - .Define("xdetf","xdet["+str_i+"]") - .Define("zdetf","zdet["+str_i+"]") - .Define("rotationf","rotation["+str_i+"]"); - - // auto report = filterDF.Display({"NHits","theta","energy"}); - // report->Print(); - //Calculate Min and Max values - auto xminf = filterDF.Min("xposf").GetValue(); - auto xmaxf = filterDF.Max("xposf").GetValue(); - auto yminf = filterDF.Min("yposf").GetValue(); - auto ymaxf = filterDF.Max("yposf").GetValue(); - auto pxminf = filterDF.Min("xmomf").GetValue(); - auto pxmaxf = filterDF.Max("xmomf").GetValue(); - auto pyminf = filterDF.Min("ymomf").GetValue(); - auto pymaxf = filterDF.Max("ymomf").GetValue(); - // Calculate means and standard deviations - xMeans[name] = filterDF.Mean("xposf").GetValue(); - yMeans[name] = filterDF.Mean("yposf").GetValue(); - xStdDevs[name] = filterDF.StdDev("xposf").GetValue(); - yStdDevs[name] = filterDF.StdDev("yposf").GetValue(); - pxMeans[name] = filterDF.Mean("xmomf").GetValue(); - pyMeans[name] = filterDF.Mean("ymomf").GetValue(); - pxStdDevs[name] = filterDF.StdDev("xmomf").GetValue(); - pyStdDevs[name] = filterDF.StdDev("ymomf").GetValue(); - + .Define("pipeRadiusf","pipeRadius["+str_i+"]"); + TString beamspotName = "Beamspot ID"+str_i+";x offset [cm]; y offset [cm]"; TString xyname = name+";x Offset [cm]; y Offset [cm]"; @@ -259,9 +160,6 @@ int acceptanceAnalysis( TString inFile = "/scratch/EIC/G4out/beamlin //Parameters of the pipe pipeRadii[name] = filterDF.Max("pipeRadiusf").GetValue(); std::cout << "Pipe ID: " << name << " Radius: " << pipeRadii[name] << " " << filterDF.Min("pipeRadiusf").GetValue() << std::endl; - pipeXPos[name] = filterDF.Max("xdetf").GetValue(); - pipeZPos[name] = filterDF.Max("zdetf").GetValue(); - pipeRotation[name] = filterDF.Max("rotationf").GetValue(); } From 96e8b507aadda7f1f474cb50be5a9ce30525a097 Mon Sep 17 00:00:00 2001 From: simonge Date: Mon, 23 Jun 2025 23:12:04 +0100 Subject: [PATCH 28/36] Define both simulations in the yml --- benchmarks/beamline/config.yml | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/benchmarks/beamline/config.yml b/benchmarks/beamline/config.yml index ef6b84e4..f92d9e7c 100644 --- a/benchmarks/beamline/config.yml +++ b/benchmarks/beamline/config.yml @@ -1,4 +1,4 @@ -sim:beamline: +sim:beamline:beamspot: extends: .det_benchmark stage: simulate script: @@ -6,6 +6,14 @@ sim:beamline: snakemake $SNAKEMAKE_FLAGS --cores 2 \ sim_output/beamline/beamlineTestlocal.edm4hep.root +sim:beamline:acceptance: + extends: .det_benchmark + stage: simulate + script: + - | + snakemake $SNAKEMAKE_FLAGS --cores 2 \ + sim_output/beamline/acceptanceTestlocal.edm4hep.root + bench:beamline: extends: .det_benchmark stage: benchmarks From 6878af8cba5902da8c27171e6d8aaaf2fdcb491e Mon Sep 17 00:00:00 2001 From: simonge Date: Tue, 24 Jun 2025 12:47:06 +0100 Subject: [PATCH 29/36] Add entry fraction plot --- benchmarks/beamline/Snakefile | 7 ++-- benchmarks/beamline/acceptanceAnalysis.C | 44 ++++++++++++++++++++---- 2 files changed, 43 insertions(+), 8 deletions(-) diff --git a/benchmarks/beamline/Snakefile b/benchmarks/beamline/Snakefile index a1cc0911..0d40f90a 100644 --- a/benchmarks/beamline/Snakefile +++ b/benchmarks/beamline/Snakefile @@ -68,9 +68,11 @@ rule beamline_acceptance_analysis: beampipe_canvas=ANALYSISDIR+"acceptance_in_beampipe_{CAMPAIGN}.png", etheta_canvas=ANALYSISDIR+"acceptance_energy_theta_{CAMPAIGN}.png", etheta_acceptance_canvas=ANALYSISDIR+"acceptance_energy_theta_acceptance_{CAMPAIGN}.png", + entrys_canvas=ANALYSISDIR+"acceptance_entrys_{CAMPAIGN}.png", shell: """ - root -l -b -q '{input.script}("{input.data}", "{output.rootfile}", "{params.xml}", "{output.beampipe_canvas}","{output.etheta_canvas}","{output.etheta_acceptance_canvas}")' + root -l -b -q '{input.script}("{input.data}", "{output.rootfile}", "{params.xml}", "{output.beampipe_canvas}","{output.etheta_canvas}","{output.etheta_acceptance_canvas}", + "{output.entrys_canvas}")' """ rule beamline: @@ -85,7 +87,8 @@ rule beamline: ANALYSISDIR+"acceptanceTestAnalysis{CAMPAIGN}.root", ANALYSISDIR+"acceptance_in_beampipe_{CAMPAIGN}.png", ANALYSISDIR+"acceptance_energy_theta_{CAMPAIGN}.png", - ANALYSISDIR+"acceptance_energy_theta_acceptance_{CAMPAIGN}.png", + ANALYSISDIR+"acceptance_energy_theta_acceptance_{CAMPAIGN}.png", + ANALYSISDIR+"acceptance_entrys_{CAMPAIGN}.png", output: directory("results/beamline/{CAMPAIGN}/") shell: diff --git a/benchmarks/beamline/acceptanceAnalysis.C b/benchmarks/beamline/acceptanceAnalysis.C index 19cb21ff..0cf6c72e 100644 --- a/benchmarks/beamline/acceptanceAnalysis.C +++ b/benchmarks/beamline/acceptanceAnalysis.C @@ -22,9 +22,10 @@ using RNode = ROOT::RDF::RNode; int acceptanceAnalysis( TString inFile = "/scratch/EIC/G4out/beamline/acceptanceTest.edm4hep.root", TString outFile = "output.root", std::string compactName = "/home/simong/EIC/epic/install/share/epic/epic_ip6_extended.xml", - TString beampipeCanvasName = "phasespace_in_beampipe.png", - TString EThetaCanvasName = "phasespace_energy_theta.png", - TString EThetaAccCanvasName= "phasespace_energy_theta_acceptance.png") { + TString beampipeCanvasName = "acceptance_in_beampipe.png", + TString EThetaCanvasName = "acceptance_energy_theta.png", + TString EThetaAccCanvasName= "acceptance_energy_theta_acceptance.png", + TString entryFractionCanvasName = "acceptance_entries.png") { //Set ROOT style gStyle->SetPadLeftMargin(0.1); // Set left margin @@ -52,6 +53,8 @@ int acceptanceAnalysis( TString inFile = "/scratch/EIC/G4out/beamlin RNode d1 = d0; RVecS colNames = d1.GetColumnNames(); + //Get total number of entries + int nEntries = d1.Count().GetValue(); //Set number of entries to process // d1 = d1.Range(0,1000); @@ -92,7 +95,10 @@ int acceptanceAnalysis( TString inFile = "/scratch/EIC/G4out/beamlin return 1; } - auto totalETheta = d1.Histo2D({"Energy vs Theta","Energy vs Theta",100,4,18,100,0,0.011},"energy","theta"); + int eBins = 100; + int thetaBins = 100; + + auto totalETheta = d1.Histo2D({"Energy vs Theta","Energy vs Theta",eBins,4,18,thetaBins,0,0.011},"energy","theta"); if(Any(colNames==readoutName)){ @@ -133,6 +139,8 @@ int acceptanceAnalysis( TString inFile = "/scratch/EIC/G4out/beamlin std::map pipeRadii; + std::map filterEntries; + std::map filterAcceptanceIntegral; //Create histograms for(int i=0; i<=7; i++){ @@ -155,7 +163,7 @@ int acceptanceAnalysis( TString inFile = "/scratch/EIC/G4out/beamlin auto extraFilterDF = filterDF.Filter(std::to_string(i+1)+"==NHits" ); TString EThetaName = "Energy vs Theta ID"+str_i+";Energy [GeV]; Theta [rad]"; TString nameETheta = name+";Energy [GeV]; Theta [rad]"; - hHistsETheta[name] = extraFilterDF.Histo2D({EThetaName,EThetaName,100,4,18,100,0,0.011},"energy","theta"); + hHistsETheta[name] = extraFilterDF.Histo2D({EThetaName,EThetaName,eBins,4,18,thetaBins,0,0.011},"energy","theta"); //Parameters of the pipe pipeRadii[name] = filterDF.Max("pipeRadiusf").GetValue(); @@ -186,9 +194,10 @@ int acceptanceAnalysis( TString inFile = "/scratch/EIC/G4out/beamlin TCanvas *cETheta = new TCanvas("energy_theta_canvas","energy_theta_canvas",3000,1600); cETheta->Divide(4,2); i=1; - for(auto [name,h] : hHistsETheta){ + for(auto [name,h] : hHistsETheta){ cETheta->cd(i++); h->Draw("colz"); + filterEntries[name] = h->GetEntries()/ nEntries; } // Canvas for energy vs theta acceptance @@ -199,17 +208,40 @@ int acceptanceAnalysis( TString inFile = "/scratch/EIC/G4out/beamlin cEThetaAcc->cd(i++); h->Divide(totalETheta.GetPtr()); h->Draw("colz"); + filterAcceptanceIntegral[name] = h->Integral()/ (eBins*thetaBins); // Normalize by the number of bins } + TH1F* hPipeEntries = CreateFittedHistogram("hNumberOfHits", + "Number of Hits per Pipe ID", + filterEntries, + {}, + "Pipe ID"); + TH1F* hPipeAcceptance = CreateFittedHistogram("hPipeAcceptance", + "Acceptance per Pipe ID", + filterAcceptanceIntegral, + {}, + "Pipe ID"); + + TCanvas *cPipeAcceptance = new TCanvas("cPipeAcceptance", "Pipe Acceptance", 1200, 400); + cPipeAcceptance->Divide(2, 1); + cPipeAcceptance->cd(1); + hPipeEntries->Draw(""); + cPipeAcceptance->cd(2); + hPipeAcceptance->Draw(""); + cPipeAcceptance->SetGrid(); + cPipeAcceptance->Update(); + // Save 2D canvases as pngs cXY->SaveAs(beampipeCanvasName); cETheta->SaveAs(EThetaCanvasName); cEThetaAcc->SaveAs(EThetaAccCanvasName); + cPipeAcceptance->SaveAs(entryFractionCanvasName); TFile *f = new TFile(outFile,"RECREATE"); cXY->Write(); cETheta->Write(); cEThetaAcc->Write(); + cPipeAcceptance->Write(); f->Close(); From 09f8681f8b255ede321b18769d74b69e95bd5c6f Mon Sep 17 00:00:00 2001 From: simonge Date: Tue, 24 Jun 2025 19:46:39 +0100 Subject: [PATCH 30/36] Make filtering more robust --- benchmarks/beamline/acceptanceAnalysis.C | 19 +++++----- benchmarks/beamline/beamlineAnalysis.C | 47 +++++++++++++++--------- 2 files changed, 39 insertions(+), 27 deletions(-) diff --git a/benchmarks/beamline/acceptanceAnalysis.C b/benchmarks/beamline/acceptanceAnalysis.C index 0cf6c72e..5f530408 100644 --- a/benchmarks/beamline/acceptanceAnalysis.C +++ b/benchmarks/beamline/acceptanceAnalysis.C @@ -138,7 +138,7 @@ int acceptanceAnalysis( TString inFile = "/scratch/EIC/G4out/beamlin std::map> hHistsETheta; - std::map pipeRadii; + std::map> pipeRadii; std::map filterEntries; std::map filterAcceptanceIntegral; @@ -148,10 +148,9 @@ int acceptanceAnalysis( TString inFile = "/scratch/EIC/G4out/beamlin std::string name = "pipeID"; std::string str_i = std::to_string(i); name += str_i; - auto filterDF = d1.Filter(std::to_string(i+1)+"<=NHits" ) - .Define("xposf","xpos["+str_i+"]") - .Define("yposf","ypos["+str_i+"]") - .Define("pipeRadiusf","pipeRadius["+str_i+"]"); + auto filterDF = d1.Define("xposf","xpos[pipeID=="+str_i+"]") + .Define("yposf","ypos[pipeID=="+str_i+"]") + .Define("pipeRadiusf","pipeRadius[pipeID=="+str_i+"]"); TString beamspotName = "Beamspot ID"+str_i+";x offset [cm]; y offset [cm]"; @@ -160,15 +159,15 @@ int acceptanceAnalysis( TString inFile = "/scratch/EIC/G4out/beamlin TString yname = name+";y Offset [cm]; y trajectory component"; hHistsxy[name] = filterDF.Histo2D({beamspotName,beamspotName,400,-maxPipeRadius,maxPipeRadius,400,-maxPipeRadius,maxPipeRadius},"xposf","yposf"); - auto extraFilterDF = filterDF.Filter(std::to_string(i+1)+"==NHits" ); + auto extraFilterDF = filterDF.Filter("All(pipeID<"+std::to_string(i+1)+")&&Any(pipeID=="+std::to_string(i)+")" ); TString EThetaName = "Energy vs Theta ID"+str_i+";Energy [GeV]; Theta [rad]"; TString nameETheta = name+";Energy [GeV]; Theta [rad]"; hHistsETheta[name] = extraFilterDF.Histo2D({EThetaName,EThetaName,eBins,4,18,thetaBins,0,0.011},"energy","theta"); //Parameters of the pipe - pipeRadii[name] = filterDF.Max("pipeRadiusf").GetValue(); - std::cout << "Pipe ID: " << name << " Radius: " << pipeRadii[name] << " " << filterDF.Min("pipeRadiusf").GetValue() << std::endl; - + pipeRadii[name] = filterDF.Max("pipeRadiusf"); + // std::cout << "Pipe ID: " << name << " Radius: " << pipeRadii[name] << " " << filterDF.Min("pipeRadiusf").GetValue() << std::endl; + } // Create histograms of the beamspot @@ -177,7 +176,7 @@ int acceptanceAnalysis( TString inFile = "/scratch/EIC/G4out/beamlin int i=1; for(auto [name,h] : hHistsxy){ // Get the pipe radius for this histogram - auto pipeRadius = pipeRadii[name]; + auto pipeRadius = pipeRadii[name].GetValue(); cXY->cd(i++); h->Draw("col"); diff --git a/benchmarks/beamline/beamlineAnalysis.C b/benchmarks/beamline/beamlineAnalysis.C index 851c4fd0..f83196cc 100644 --- a/benchmarks/beamline/beamlineAnalysis.C +++ b/benchmarks/beamline/beamlineAnalysis.C @@ -59,7 +59,7 @@ int beamlineAnalysis( TString inFile = "/scratch/EIC/G4out/beamline/b //Set number of entries to process // d1 = d1.Range(0,1000); int nEntries = d1.Count().GetValue(); - float acceptableLoss = 0.99; // Set the acceptable loss percentage to 1% + float acceptableLoss = 0.999; // Set the acceptable loss percentage to 0.1% float acceptableEntries = nEntries * acceptableLoss; //Get the collection @@ -171,28 +171,39 @@ int beamlineAnalysis( TString inFile = "/scratch/EIC/G4out/beamline/b std::map pipeZPos; std::map pipeRotation; - auto xmin = d1.Min("xpos").GetValue(); - auto xmax = d1.Max("xpos").GetValue(); - auto ymin = d1.Min("ypos").GetValue(); - auto ymax = d1.Max("ypos").GetValue(); - auto pxmin = d1.Min("xmom").GetValue(); - auto pxmax = d1.Max("xmom").GetValue(); - auto pymin = d1.Min("ymom").GetValue(); - auto pymax = d1.Max("ymom").GetValue(); + // Queue up all actions + auto xmin_ptr = d1.Min("xpos"); + auto xmax_ptr = d1.Max("xpos"); + auto ymin_ptr = d1.Min("ypos"); + auto ymax_ptr = d1.Max("ypos"); + auto pxmin_ptr = d1.Min("xmom"); + auto pxmax_ptr = d1.Max("xmom"); + auto pymin_ptr = d1.Min("ymom"); + auto pymax_ptr = d1.Max("ymom"); + + // Now trigger the event loop (only once) + auto xmin = xmin_ptr.GetValue(); + auto xmax = xmax_ptr.GetValue(); + auto ymin = ymin_ptr.GetValue(); + auto ymax = ymax_ptr.GetValue(); + auto pxmin = pxmin_ptr.GetValue(); + auto pxmax = pxmax_ptr.GetValue(); + auto pymin = pymin_ptr.GetValue(); + auto pymax = pymax_ptr.GetValue(); //Create histograms for(int i=0; i<=7; i++){ std::string name = "pipeID"; name += std::to_string(i); - auto filterDF = d1.Define("xposf","xpos["+std::to_string(i)+"]") - .Define("yposf","ypos["+std::to_string(i)+"]") - .Define("xmomf","xmom["+std::to_string(i)+"]") - .Define("ymomf","ymom["+std::to_string(i)+"]") - .Define("pipeRadiusf","pipeRadius["+std::to_string(i)+"]") - .Define("xdetf","xdet["+std::to_string(i)+"]") - .Define("zdetf","zdet["+std::to_string(i)+"]") - .Define("rotationf","rotation["+std::to_string(i)+"]"); + auto filterDF = d1.Define("xposf","xpos[pipeID=="+std::to_string(i)+"]") + .Define("yposf","ypos[pipeID=="+std::to_string(i)+"]") + .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("xdetf","xdet[pipeID=="+std::to_string(i)+"]") + .Define("zdetf","zdet[pipeID=="+std::to_string(i)+"]") + .Define("rotationf","rotation[pipeID=="+std::to_string(i)+"]"); //Calculate Min and Max values @@ -292,6 +303,8 @@ int beamlineAnalysis( TString inFile = "/scratch/EIC/G4out/beamline/b } + + // Create histograms of the beamspot TCanvas *cXY = new TCanvas("beamspot_canvas","beamspot_canvas",3000,1600); cXY->Divide(4,2); From 40e1fa9a881ceb3513f9d376113aab052c4d2210 Mon Sep 17 00:00:00 2001 From: simonge Date: Tue, 24 Jun 2025 20:49:20 +0100 Subject: [PATCH 31/36] Change entry limit warning --- benchmarks/beamline/beamlineAnalysis.C | 2 +- benchmarks/beamline/beamlineGPS.mac | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/benchmarks/beamline/beamlineAnalysis.C b/benchmarks/beamline/beamlineAnalysis.C index f83196cc..42525598 100644 --- a/benchmarks/beamline/beamlineAnalysis.C +++ b/benchmarks/beamline/beamlineAnalysis.C @@ -314,7 +314,7 @@ int beamlineAnalysis( TString inFile = "/scratch/EIC/G4out/beamline/b // Check histogram entries if(h->GetEntries() < acceptableEntries){ std::cout << "Warning: Only " << h->GetEntries()/nEntries << " of particles contributing to histogram " << name - << " , which is below the accepted threshold of " << acceptableEntries << std::endl; + << " , which is below the accepted threshold of " << acceptableEntries/nEntries << std::endl; pass = 1; } diff --git a/benchmarks/beamline/beamlineGPS.mac b/benchmarks/beamline/beamlineGPS.mac index d3cdc21b..54267451 100644 --- a/benchmarks/beamline/beamlineGPS.mac +++ b/benchmarks/beamline/beamlineGPS.mac @@ -8,4 +8,4 @@ /gps/energy 18 GeV /gps/particle e- -/run/beamOn 10000 \ No newline at end of file +/run/beamOn 100000 \ No newline at end of file From 158c8efb220eacefb6c407d478ed117c57db6d65 Mon Sep 17 00:00:00 2001 From: simonge Date: Wed, 25 Jun 2025 13:21:39 +0100 Subject: [PATCH 32/36] Add reconstruction training based on beampipe exit --- benchmarks/beamline/ProcessData.py | 15 ++ benchmarks/beamline/RegressionModel.py | 128 +++++++++++++ benchmarks/beamline/Snakefile | 12 ++ benchmarks/beamline/SteeringRegression.py | 55 ++++++ benchmarks/beamline/TestModel.py | 211 ++++++++++++++++++++++ benchmarks/beamline/processData.C | 37 ++++ 6 files changed, 458 insertions(+) create mode 100644 benchmarks/beamline/ProcessData.py create mode 100644 benchmarks/beamline/RegressionModel.py create mode 100644 benchmarks/beamline/SteeringRegression.py create mode 100644 benchmarks/beamline/TestModel.py create mode 100644 benchmarks/beamline/processData.C diff --git a/benchmarks/beamline/ProcessData.py b/benchmarks/beamline/ProcessData.py new file mode 100644 index 00000000..7eb6ad1a --- /dev/null +++ b/benchmarks/beamline/ProcessData.py @@ -0,0 +1,15 @@ +import uproot +import awkward as ak + +def create_arrays(dataFiles,beamEnergy=18): + + # List of branches to load + branches = ["features","targets"] + + # Load data from concatenated list of files + data = uproot.concatenate([f"{file}:events" for file in dataFiles], branches, library="ak") + + input_data = data["features"] + target_data = data["targets"]/beamEnergy + + return input_data, target_data \ No newline at end of file diff --git a/benchmarks/beamline/RegressionModel.py b/benchmarks/beamline/RegressionModel.py new file mode 100644 index 00000000..ecef920a --- /dev/null +++ b/benchmarks/beamline/RegressionModel.py @@ -0,0 +1,128 @@ +import torch +import torch.nn as nn +import torch.optim as optim +import numpy as np + +class ProjectToX0Plane(nn.Module): + def forward(self, x): + # x shape: (batch, 6) -> [x, y, z, px, py, pz] + x0 = x[:, 0] + y0 = x[:, 1] + z0 = x[:, 2] + px = x[:, 3] + py = x[:, 4] + pz = x[:, 5] + + # Avoid division by zero for px + eps = 1e-8 + px_safe = torch.where(px.abs() < eps, eps * torch.sign(px) + eps, px) + t = -x0 / px_safe + + y_proj = y0 + py * t + z_proj = z0 + pz * t + + # Output: [y_proj, z_proj, px, pz] + return torch.stack([y_proj, z_proj, px, pz], dim=1) + +class RegressionModel(nn.Module): + def __init__(self): + super(RegressionModel, self).__init__() + self.fc1 = nn.Linear(4, 512) + self.fc2 = nn.Linear(512, 64) + self.fc4 = nn.Linear(64, 3) + self.input_mean = torch.tensor([0.0, 0.0, 0.0, 0.0]) + self.input_std = torch.tensor([1.0, 1.0, 1.0, 1.0]) + # self.input_covariance = torch.tensor([[1.0, 0.0, 0.0, 0.0], + # [0.0, 1.0, 0.0, 0.0], + # [0.0, 0.0, 1.0, 0.0], + # [0.0, 0.0, 0.0, 1.0]]) + self.output_mean = torch.tensor([0.0, 0.0, 0.0]) + self.output_std = torch.tensor([1.0, 1.0, 1.0]) + # self.output_correlation = torch.tensor([[1.0, 0.0, 0.0], + # [0.0, 1.0, 0.0], + # [0.0, 0.0, 1.0]]) + + def forward(self, x): + x = ProjectToX0Plane()(x) + x = (x-self.input_mean)/self.input_std + x = torch.tanh(self.fc1(x)) + x = torch.tanh(self.fc2(x)) + x = self.fc4(x) + x = x*self.output_std + self.output_mean + return x + + def adapt(self, input_data, output_data): + in_mean = input_data.mean(axis=0) + in_std = input_data.std (axis=0) + self.input_mean = torch.tensor(in_mean) + self.input_std = torch.tensor(in_std) + + # Calculate the correlation matrix of the input data + # input_normalized = (input_data-in_mean)/in_std + # input_correlation = np.corrcoef(input_normalized, rowvar=False) + # Invert the correlation matrix and convert into float tensor + # self.input_covariance = torch.tensor(np.linalg.inv(input_correlation).astype(np.float32)) + + self.output_mean = torch.tensor(output_data.mean(axis=0)) + self.output_std = torch.tensor(output_data.std (axis=0)) + +def makeModel(): + # Create the model + model = RegressionModel() + # Define the optimizer + optimizer = optim.Adam(model.parameters(), lr=0.0001) + # Define the loss function + criterion = nn.MSELoss() + + return model, optimizer, criterion + +def trainModel(epochs, train_loader, val_loader): + + # device = torch.device("cuda" if torch.cuda.is_available() else "cpu") + # print(f"Using device: {device}") + + model, optimizer, criterion = makeModel() + # model.to(device) + + # Verify that the model parameters are on the GPU + # for name, param in model.named_parameters(): + # print(f"{name} is on {param.device}") + + # Adapt the model using the training data from the training loader + # Project inputs to X0 plane before adaptation + inputs = train_loader.dataset.tensors[0] + projected_inputs = ProjectToX0Plane()(inputs) + model.adapt(projected_inputs.detach().numpy(), train_loader.dataset.tensors[1].detach().numpy()) + + for epoch in range(epochs): + model.train() + running_loss = 0.0 + for inputs, targets in train_loader: + # inputs, targets = inputs.to(device), targets.to(device) + optimizer.zero_grad() + outputs = model(inputs) + loss = criterion(outputs, targets) + loss.backward() + optimizer.step() + running_loss += loss.item() * inputs.size(0) + + epoch_loss = running_loss / len(train_loader.dataset) + # print(f"Epoch [{epoch+1}/{epochs}], Loss: {epoch_loss:.4f}") + + + # Validation step + model.eval() + val_loss = 0.0 + with torch.no_grad(): + for val_inputs, val_targets in val_loader: + # val_inputs, val_targets = val_inputs.to(device), val_targets.to(device) + val_outputs = model(val_inputs) + val_loss += criterion(val_outputs, val_targets).item() * val_inputs.size(0) + # val_outputs = model(val_input) + # val_loss = criterion(val_outputs, val_target) + + val_loss /= len(val_loader.dataset) + + print(f"Epoch [{epoch+1}/{epochs}], Loss: {loss}, Val Loss: {val_loss}") + + return model \ No newline at end of file diff --git a/benchmarks/beamline/Snakefile b/benchmarks/beamline/Snakefile index 0d40f90a..4a0ccfb2 100644 --- a/benchmarks/beamline/Snakefile +++ b/benchmarks/beamline/Snakefile @@ -75,6 +75,18 @@ rule beamline_acceptance_analysis: "{output.entrys_canvas}")' """ +# Trains a regression model to predict the TaggerTrackerTargetTensor from the TaggerTrackerFeatureTensor. +rule beamline_steering_reconstruction_training: + input: + script="SteeringRegression.py", + data=SIMOUTDIR+"acceptanceTest{CAMPAIGN}.edm4hep.root", + output: + onnxfile=ANALYSISDIR+"BeamlineSteeringReconstruction.onnx", + shell: + """ + python {input.script} --dataFiles {input.data} --outModelFile {output.onnxfile} + """ + rule beamline: input: ANALYSISDIR+"beamlineTestAnalysis{CAMPAIGN}.root", diff --git a/benchmarks/beamline/SteeringRegression.py b/benchmarks/beamline/SteeringRegression.py new file mode 100644 index 00000000..50d30004 --- /dev/null +++ b/benchmarks/beamline/SteeringRegression.py @@ -0,0 +1,55 @@ +import torch +import argparse +from ProcessData import create_arrays +from torch.utils.data import DataLoader, TensorDataset +from RegressionModel import makeModel, trainModel + +# Parse arguments +parser = argparse.ArgumentParser(description='Train a regression model for the Tagger.') +parser.add_argument('--dataFiles', type=str, nargs='+', help='Path to the data files') +parser.add_argument('--outModelFile', type=str, default="regression_model.onnx", help='Output file for the trained model') +parser.add_argument('--batchSize', type=int, default=4096, help='Batch size for training') +parser.add_argument('--epochs', type=int, default=100, help='Number of epochs for training') +args = parser.parse_args() + +input_data, target_data = create_arrays(args.dataFiles) + +# print(f"Input data shape: {input_data.shape}") +# print(f"Target data shape: {target_data.shape}") + +torch_input_data = torch.tensor(input_data) +torch_target_data = torch.tensor(target_data) + +print(f"Input data shape: {torch_input_data.shape}") +print(f"Target data shape: {torch_target_data.shape}") + +# Split data into training and validation sets +validation_fraction = 0.25 +split_index = int(len(torch_input_data) * (1 - validation_fraction)) + +val_input_data = torch_input_data[split_index:] +val_target_data = torch_target_data[split_index:] +train_input_data = torch_input_data[:split_index] +train_target_data = torch_target_data[:split_index] + +# Create TensorDatasets +train_dataset = TensorDataset(train_input_data, train_target_data) +val_dataset = TensorDataset(val_input_data, val_target_data) + +# Create DataLoaders +train_loader = DataLoader(train_dataset, batch_size=args.batchSize, shuffle=True ) +val_loader = DataLoader(val_dataset, batch_size=args.batchSize, shuffle=False) + +print(f"Training data: {len(train_input_data)} samples") + +model = trainModel(args.epochs, train_loader, val_loader) + +# Save the trained model to ONNX format + +dummy_input = torch_input_data[0].unsqueeze(0) # Create a dummy input for the model + +torch.onnx.export(model, dummy_input, args.outModelFile, + input_names=['input'], output_names=['output'], + dynamic_axes={'input': {0: 'batch_size'}, 'output': {0: 'batch_size'}}) + +print(f"Model has been saved to {args.outModelFile}") \ No newline at end of file diff --git a/benchmarks/beamline/TestModel.py b/benchmarks/beamline/TestModel.py new file mode 100644 index 00000000..dfd95931 --- /dev/null +++ b/benchmarks/beamline/TestModel.py @@ -0,0 +1,211 @@ +import onnxruntime as ort +import argparse +import numpy as np +from ProcessData import create_arrays +import matplotlib.pyplot as plt +from matplotlib.colors import LogNorm +from scipy.stats import norm +from scipy.optimize import curve_fit + +# Parse arguments +parser = argparse.ArgumentParser(description='Train a regression model for the Tagger.') +parser.add_argument('--modelFile', type=str, default="regression_model.onnx", help='Path to the ONNX model file') +parser.add_argument('--dataFiles', type=str, nargs='+', help='Path to the data files') +parser.add_argument('--outDir', type=str, default=".", help='Output directory') +args = parser.parse_args() +modelFile = args.modelFile +dataFiles = args.dataFiles +outDir = args.outDir +outGraphFile = outDir + "/output_vs_target.png" +outGraphFile2 = outDir + "/output_vs_target2.png" +outGraphFile3 = outDir + "/output_vs_target3.png" +outGraphFile4 = outDir + "/output_vs_target4.png" + +input_data, target_data = create_arrays(dataFiles) + +target_data = np.array(target_data) + +# Load the ONNX model +session = ort.InferenceSession(modelFile) + +# Run the model on the input data +input_name = session.get_inputs()[0].name +output_name = session.get_outputs()[0].name +input_data = np.array(input_data,dtype=np.float32) +output = session.run([output_name], {input_name: input_data}) +output = np.array(output[0]) + +out_theta = np.arctan2(np.sqrt(output[:,0]**2 + output[:,1]**2),output[:,2]) +out_phi = np.arctan2(output[:,1],output[:,0]) +out_mag = np.sqrt(output[:,0]**2 + output[:,1]**2 + output[:,2]**2) +in_theta = np.arctan2(np.sqrt(target_data[:,0]**2 + target_data[:,1]**2),target_data[:,2]) +in_phi = np.arctan2(target_data[:,1],target_data[:,0]) +in_mag = np.sqrt(target_data[:,0]**2 + target_data[:,1]**2 + target_data[:,2]**2) + + +thetadiff = out_theta - in_theta +phidiff = out_phi - in_phi +# Move phidiff to within -pi/2 and pi/2 +phidiff = (phidiff + np.pi) % (2 * np.pi) - np.pi +magdiff = out_mag - in_mag + +diff = (target_data - output)/target_data +diffrange = [[-5,5],[-5,5],[-0.5,0.5]] +datarange = [[-0.02,0.02],[-0.02,0.02],[-1,0]] + +# Use the 'seismic' colormap +cmap = plt.get_cmap('seismic') + +# Creates histograms to compare the target and output data +fig, axs = plt.subplots(3, 3, figsize=(12, 12)) +for i in range(3): + # 2D histograms showing trends in the data + axs[0,i].hist2d(target_data[:,i], output[:,i], bins=100, range=[datarange[i],datarange[i]], cmap="seismic", norm=LogNorm(), label="Output vs Target") + axs[0,i].set_xlabel(f"Variable {i} Target") + axs[0,i].set_ylabel(f"Variable {i} Output") + + axs[1,i].hist(diff[:,i], bins=100, alpha=0.5, range=diffrange[i], label="Difference") + axs[1,i].set_xlabel(f"Variable {i} Difference") + axs[1,i].set_ylabel("Counts") + + axs[2,i].hist2d(target_data[:,i], diff[:,i], bins=100, range=[datarange[i],diffrange[i]], cmap="seismic", norm=LogNorm(), label="Difference vs Target") + axs[2,i].set_xlabel(f"Variable {i} Target") + axs[2,i].set_ylabel(f"Variable {i} Difference") + +plt.show() +plt.savefig(outGraphFile) + +# Creates histograms to compare theta, phi and mag target and output data +fig2, axs2 = plt.subplots(3, 3, figsize=(12, 12)) + +thetarange = [np.pi-0.01,np.pi] +phirange = [-np.pi,np.pi] +magrange = [0,1] + +thetadiffrange = [-0.01,0.01] +phidiffrange = [-np.pi,np.pi] +magdiffrange = [-0.1,0.1] + +# 2D histograms showing trends in the data +axs2[0,0].hist2d(out_theta, in_theta, bins=100, range=[thetarange,thetarange], cmap="seismic", norm=LogNorm(), label="Output vs Target") +axs2[0,0].set_xlabel("Theta Target") +axs2[0,0].set_ylabel("Theta Output") + +axs2[0,1].hist2d(out_phi, in_phi, bins=100, range=[phirange,phirange], cmap="seismic", label="Output vs Target") +axs2[0,1].set_xlabel("Phi Target") +axs2[0,1].set_ylabel("Phi Output") + +axs2[0,2].hist2d(out_mag, in_mag, bins=100, range=[magrange,magrange], cmap="seismic", norm=LogNorm(), label="Output vs Target") +axs2[0,2].set_xlabel("Mag Target") +axs2[0,2].set_ylabel("Mag Output") + +axs2[1,0].hist(thetadiff, bins=100, alpha=0.5, range=thetadiffrange, label="Difference") +axs2[1,0].set_xlabel("Theta Difference") +axs2[1,0].set_ylabel("Counts") + +axs2[1,1].hist(phidiff, bins=100, alpha=0.5, range=phidiffrange, label="Difference") +axs2[1,1].set_xlabel("Phi Difference") +axs2[1,1].set_ylabel("Counts") + +axs2[1,2].hist(magdiff, bins=100, alpha=0.5, range=magdiffrange, label="Difference") +axs2[1,2].set_xlabel("Mag Difference") +axs2[1,2].set_ylabel("Counts") + +axs2[2,0].hist2d(in_theta, thetadiff, bins=100, range=[thetarange,thetadiffrange], cmap="seismic", norm=LogNorm(), label="Difference vs Target") +axs2[2,0].set_xlabel("Theta Target") +axs2[2,0].set_ylabel("Theta Difference") + +axs2[2,1].hist2d(in_phi, phidiff, bins=100, range=[phirange,phidiffrange], cmap="seismic", label="Difference vs Target") +axs2[2,1].set_xlabel("Phi Target") +axs2[2,1].set_ylabel("Phi Difference") + +axs2[2,2].hist2d(in_mag, magdiff, bins=100, range=[magrange,magdiffrange], cmap="seismic", norm=LogNorm(), label="Difference vs Target") +axs2[2,2].set_xlabel("Mag Target") +axs2[2,2].set_ylabel("Mag Difference") + +plt.show() +plt.savefig(outGraphFile2) + +# Create histograms where the theta value has been cut at less than 3.14 +fig3, axs3 = plt.subplots(3, 3, figsize=(12, 12)) + +out_theta_cut = out_theta[out_theta < 3.14] +in_theta_cut = in_theta[out_theta < 3.14] +thetadiff_cut = thetadiff[out_theta < 3.14] + +out_phi_cut = out_phi[out_theta < 3.14] +in_phi_cut = in_phi[out_theta < 3.14] +phidiff_cut = phidiff[out_theta < 3.14] + +out_mag_cut = out_mag[out_theta < 3.14] +in_mag_cut = in_mag[out_theta < 3.14] +magdiff_cut = magdiff[out_theta < 3.14] + +axs3[0,0].hist2d(out_theta_cut, in_theta_cut, bins=100, range=[thetarange,thetarange], cmap="seismic", norm=LogNorm(), label="Output vs Target") +axs3[0,0].set_xlabel("Theta Target") +axs3[0,0].set_ylabel("Theta Output") + +axs3[0,1].hist2d(out_phi_cut, in_phi_cut, bins=100, range=[phirange,phirange], cmap="seismic", label="Output vs Target") +axs3[0,1].set_xlabel("Phi Target") +axs3[0,1].set_ylabel("Phi Output") + +axs3[0,2].hist2d(out_mag_cut, in_mag_cut, bins=100, range=[magrange,magrange], cmap="seismic", norm=LogNorm(), label="Output vs Target") +axs3[0,2].set_xlabel("Mag Target") +axs3[0,2].set_ylabel("Mag Output") + +axs3[1,0].hist(thetadiff_cut, bins=100, alpha=0.5, range=thetadiffrange, label="Difference") +axs3[1,0].set_xlabel("Theta Difference") +axs3[1,0].set_ylabel("Counts") + +axs3[1,1].hist(phidiff_cut, bins=100, alpha=0.5, range=phidiffrange, label="Difference") +axs3[1,1].set_xlabel("Phi Difference") +axs3[1,1].set_ylabel("Counts") + +axs3[1,2].hist(magdiff_cut, bins=100, alpha=0.5, range=magdiffrange, label="Difference") +axs3[1,2].set_xlabel("Mag Difference") +axs3[1,2].set_ylabel("Counts") + +axs3[2,0].hist2d(in_theta_cut, thetadiff_cut, bins=100, range=[thetarange,thetadiffrange], cmap="seismic", norm=LogNorm(), label="Difference vs Target") +axs3[2,0].set_xlabel("Theta Target") +axs3[2,0].set_ylabel("Theta Difference") + +axs3[2,1].hist2d(in_phi_cut, phidiff_cut, bins=100, range=[phirange,phidiffrange], cmap="seismic", label="Difference vs Target") +axs3[2,1].set_xlabel("Phi Target") +axs3[2,1].set_ylabel("Phi Difference") + +axs3[2,2].hist2d(in_mag_cut, magdiff_cut, bins=100, range=[magrange,magdiffrange], cmap="seismic", norm=LogNorm(), label="Difference vs Target") +axs3[2,2].set_xlabel("Mag Target") +axs3[2,2].set_ylabel("Mag Difference") + +plt.show() +plt.savefig(outGraphFile3) + +# Create plots where a Gaussian has been fitted to the data +# Function to fit a Gaussian and plot the results +def plot_gaussian_fit(ax, data, range, xlabel, ylabel): + def gaussian(x, mu, sigma, amplitude): + return amplitude * np.exp(-0.5 * ((x - mu) / sigma) ** 2) + + hist, bin_edges = np.histogram(data, bins=100, range=range, density=True) + bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2 + + popt, _ = curve_fit(gaussian, bin_centers, hist, p0=[0, 0.01, 1]) + mu, sigma, amplitude = popt + print(f"mu={mu}, sigma={sigma}, amplitude={amplitude}") + + x = np.linspace(range[0], range[1], 100) + ax.plot(x, gaussian(x, *popt), 'k', linewidth=2) + ax.set_xlabel(xlabel) + ax.set_ylabel(ylabel) + ax.hist(data, bins=100, alpha=0.5, range=range, edgecolor='black', density=True) + ax.legend([f'Fit: $\mu$={mu:.5f}, $\sigma$={sigma:.5f}']) + +# Create histograms with Gaussian fits +fig4, axs4 = plt.subplots(3, 1, figsize=(8, 12)) + +plot_gaussian_fit(axs4[0], thetadiff, thetadiffrange, "Theta Difference", "Density") +plot_gaussian_fit(axs4[1], phidiff_cut, phidiffrange, "Phi Difference", "Density") +plot_gaussian_fit(axs4[2], magdiff, magdiffrange, "Mag Difference", "Density") + +plt.show() +plt.savefig(outGraphFile4) \ No newline at end of file diff --git a/benchmarks/beamline/processData.C b/benchmarks/beamline/processData.C new file mode 100644 index 00000000..e2262b31 --- /dev/null +++ b/benchmarks/beamline/processData.C @@ -0,0 +1,37 @@ +#include "ROOT/RDataFrame.hxx" +#include "TFile.h" +#include "TTree.h" +#include "edm4hep/MCParticle.h" +#include "edm4hep/SimTrackerHit.h" +#include + +void processData(const TString inputFile="/home/simong/EIC/detector_benchmarks_anl/sim_output/beamline/acceptanceTestcurrent.edm4hep.root", const TString outputFile="test.root", const int desired_cellID = 66757) { + // Define the branches to read + std::vector branches = { + "BackwardsBeamlineHits.cellID", + "BackwardsBeamlineHits.position.x", + "BackwardsBeamlineHits.position.y", + "BackwardsBeamlineHits.position.z", + "BackwardsBeamlineHits.momentum.x", + "BackwardsBeamlineHits.momentum.y", + "BackwardsBeamlineHits.momentum.z", + "MCParticles.momentum.x", + "MCParticles.momentum.y", + "MCParticles.momentum.z" + }; + + // Create a ROOT DataFrame to read the input files + ROOT::RDataFrame df("events", inputFile); + + //Filter on events with only a single MCParticle + auto filterDF = df.Filter("MCParticles.size()==1") + .Define("desiredHits", "BackwardsBeamlineHits[BackwardsBeamlineHits.cellID=="+std::to_string(desired_cellID)+"]") + .Filter("desiredHits.size()==1") + .Define("features", "std::array{static_cast(desiredHits[0].position.x), static_cast(desiredHits[0].position.y), static_cast(desiredHits[0].position.z), static_cast(desiredHits[0].momentum.x), static_cast(desiredHits[0].momentum.y), static_cast(desiredHits[0].momentum.z)}") + .Define("targets", "std::array{static_cast(MCParticles[0].momentum.x), static_cast(MCParticles[0].momentum.y), static_cast(MCParticles[0].momentum.z)}"); + + // Save the filtered data to a new ROOT file + filterDF.Snapshot("events", outputFile, {"features","targets"}); + + std::cout << "Filtered data saved to " << outputFile << std::endl; +} \ No newline at end of file From d937fdb0790640401847ca0203dc7c088e680efe Mon Sep 17 00:00:00 2001 From: simonge Date: Wed, 25 Jun 2025 15:14:37 +0100 Subject: [PATCH 33/36] Update model and atempt to setup snakemake --- benchmarks/beamline/RegressionModel.py | 85 +++++++++++++++----------- benchmarks/beamline/Snakefile | 35 +++++++++-- 2 files changed, 80 insertions(+), 40 deletions(-) diff --git a/benchmarks/beamline/RegressionModel.py b/benchmarks/beamline/RegressionModel.py index ecef920a..02ce775f 100644 --- a/benchmarks/beamline/RegressionModel.py +++ b/benchmarks/beamline/RegressionModel.py @@ -27,44 +27,59 @@ def forward(self, x): class RegressionModel(nn.Module): def __init__(self): super(RegressionModel, self).__init__() + self.project_to_x0 = ProjectToX0Plane() self.fc1 = nn.Linear(4, 512) self.fc2 = nn.Linear(512, 64) self.fc4 = nn.Linear(64, 3) - self.input_mean = torch.tensor([0.0, 0.0, 0.0, 0.0]) - self.input_std = torch.tensor([1.0, 1.0, 1.0, 1.0]) - # self.input_covariance = torch.tensor([[1.0, 0.0, 0.0, 0.0], - # [0.0, 1.0, 0.0, 0.0], - # [0.0, 0.0, 1.0, 0.0], - # [0.0, 0.0, 0.0, 1.0]]) - self.output_mean = torch.tensor([0.0, 0.0, 0.0]) - self.output_std = torch.tensor([1.0, 1.0, 1.0]) - # self.output_correlation = torch.tensor([[1.0, 0.0, 0.0], - # [0.0, 1.0, 0.0], - # [0.0, 0.0, 1.0]]) + + # Normalization parameters + self.input_mean = nn.Parameter(torch.zeros(4), requires_grad=False) + self.input_std = nn.Parameter(torch.ones(4), requires_grad=False) + self.output_mean = nn.Parameter(torch.zeros(3), requires_grad=False) + self.output_std = nn.Parameter(torch.ones(3), requires_grad=False) def forward(self, x): - x = ProjectToX0Plane()(x) - x = (x-self.input_mean)/self.input_std + # Apply projection and normalization + x = self.project_to_x0(x) + x = (x - self.input_mean) / self.input_std + + # Pass through the fully connected layers + x = self._core_forward(x) + + # Denormalize outputs + x = x * self.output_std + self.output_mean + return x + + def _core_forward(self, x): + # Core fully connected layers x = torch.tanh(self.fc1(x)) x = torch.tanh(self.fc2(x)) x = self.fc4(x) - x = x*self.output_std + self.output_mean return x def adapt(self, input_data, output_data): - in_mean = input_data.mean(axis=0) - in_std = input_data.std (axis=0) - self.input_mean = torch.tensor(in_mean) - self.input_std = torch.tensor(in_std) + # Compute normalization parameters from training data + self.input_mean.data = torch.tensor(input_data.mean(axis=0), dtype=torch.float32) + self.input_std.data = torch.tensor(input_data.std(axis=0), dtype=torch.float32) + self.output_mean.data = torch.tensor(output_data.mean(axis=0), dtype=torch.float32) + self.output_std.data = torch.tensor(output_data.std(axis=0), dtype=torch.float32) + +def preprocess_data(model, data_loader): + inputs = data_loader.dataset.tensors[0] + targets = data_loader.dataset.tensors[1] - # Calculate the correlation matrix of the input data - # input_normalized = (input_data-in_mean)/in_std - # input_correlation = np.corrcoef(input_normalized, rowvar=False) - # Invert the correlation matrix and convert into float tensor - # self.input_covariance = torch.tensor(np.linalg.inv(input_correlation).astype(np.float32)) + # Apply projection + projected_inputs = ProjectToX0Plane()(inputs) + + # Compute normalization parameters + model.adapt(projected_inputs.detach().numpy(), targets.detach().numpy()) + + # Normalize inputs and targets + normalized_inputs = (projected_inputs - model.input_mean) / model.input_std + normalized_targets = (targets - model.output_mean) / model.output_std - self.output_mean = torch.tensor(output_data.mean(axis=0)) - self.output_std = torch.tensor(output_data.std (axis=0)) + # Replace the dataset with preprocessed data + data_loader.dataset.tensors = (normalized_inputs, normalized_targets) def makeModel(): # Create the model @@ -78,21 +93,19 @@ def makeModel(): def trainModel(epochs, train_loader, val_loader): - # device = torch.device("cuda" if torch.cuda.is_available() else "cpu") - # print(f"Using device: {device}") + device = torch.device("cuda" if torch.cuda.is_available() else "cpu") + print(f"Using device: {device}") model, optimizer, criterion = makeModel() - # model.to(device) + model.to(device) # Verify that the model parameters are on the GPU # for name, param in model.named_parameters(): # print(f"{name} is on {param.device}") - - # Adapt the model using the training data from the training loader - # Project inputs to X0 plane before adaptation - inputs = train_loader.dataset.tensors[0] - projected_inputs = ProjectToX0Plane()(inputs) - model.adapt(projected_inputs.detach().numpy(), train_loader.dataset.tensors[1].detach().numpy()) + + # Preprocess training and validation data + preprocess_data(model, train_loader) + preprocess_data(model, val_loader) for epoch in range(epochs): model.train() @@ -100,7 +113,7 @@ def trainModel(epochs, train_loader, val_loader): for inputs, targets in train_loader: # inputs, targets = inputs.to(device), targets.to(device) optimizer.zero_grad() - outputs = model(inputs) + outputs = model._core_forward(inputs) loss = criterion(outputs, targets) loss.backward() optimizer.step() @@ -116,7 +129,7 @@ def trainModel(epochs, train_loader, val_loader): with torch.no_grad(): for val_inputs, val_targets in val_loader: # val_inputs, val_targets = val_inputs.to(device), val_targets.to(device) - val_outputs = model(val_inputs) + val_outputs = model._core_forward(val_inputs) val_loss += criterion(val_outputs, val_targets).item() * val_inputs.size(0) # val_outputs = model(val_input) # val_loss = criterion(val_outputs, val_target) diff --git a/benchmarks/beamline/Snakefile b/benchmarks/beamline/Snakefile index 4a0ccfb2..240a181d 100644 --- a/benchmarks/beamline/Snakefile +++ b/benchmarks/beamline/Snakefile @@ -75,18 +75,45 @@ rule beamline_acceptance_analysis: "{output.entrys_canvas}")' """ -# Trains a regression model to predict the TaggerTrackerTargetTensor from the TaggerTrackerFeatureTensor. -rule beamline_steering_reconstruction_training: +# Processes the simulation output data for training +rule beamline_steering_reconstruction_preperation: input: - script="SteeringRegression.py", + script=workflow.source_path("processData.C"), data=SIMOUTDIR+"acceptanceTest{CAMPAIGN}.edm4hep.root", output: - onnxfile=ANALYSISDIR+"BeamlineSteeringReconstruction.onnx", + onnxfile=ANALYSISDIR+"BeamlineSteeringReconstructionData{CAMPAIGN}.root", + shell: + """ + root -l -b -q '{input.script}("{input.data}", "{output.rootfile}")' + """ + +# Trains a regression model to predict the MCParticle momentum from the B2eR beamline virtual tracker hits. +rule beamline_steering_reconstruction_training: + input: + script=workflow.source_path("SteeringRegression.py"), + scriptdata=workflow.source_path("ProcessData.py"), + scriptmodel=workflow.source_path("RegressionModel.py"), + data=ANALYSISDIR+"BeamlineSteeringReconstructionData{CAMPAIGN}.root", + output: + onnxfile=ANALYSISDIR+"BeamlineSteeringReconstruction{CAMPAIGN}.onnx", shell: """ python {input.script} --dataFiles {input.data} --outModelFile {output.onnxfile} """ +# Trains a regression model to predict the MCParticle momentum from the B2eR beamline virtual tracker hits. +rule beamline_steering_reconstruction_training: + input: + script=workflow.source_path("TestModel.py"), + data=ANALYSISDIR+"BeamlineSteeringReconstructionData{CAMPAIGN}.root", + model=ANALYSISDIR+"BeamlineSteeringReconstruction{CAMPAIGN}.onnx", + output: + directory(ANALYSISDIR+"/NN_Test_{CAMPAIGN}/") + shell: + """ + python {input.script} --dataFiles {input.data} --modelFile {output.onnxfile} --outDir {output} + """ + rule beamline: input: ANALYSISDIR+"beamlineTestAnalysis{CAMPAIGN}.root", From e5a4ab8625e266e444a6887a9a7b1bcd62e74a1f Mon Sep 17 00:00:00 2001 From: simonge Date: Wed, 25 Jun 2025 17:31:30 +0100 Subject: [PATCH 34/36] Fix snakemane rule and silence npsim info --- benchmarks/beamline/Snakefile | 8 +++++--- benchmarks/beamline/acceptanceGPS.mac | 2 +- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/benchmarks/beamline/Snakefile b/benchmarks/beamline/Snakefile index 240a181d..32338c9a 100644 --- a/benchmarks/beamline/Snakefile +++ b/benchmarks/beamline/Snakefile @@ -13,6 +13,7 @@ rule beamline_steering_sim: --enableG4GPS \ --macroFile {input.macro} \ --compactFile $DETECTOR_PATH/epic_ip6_extended.xml \ + --printLevel WARNING \ --outputFile {output} \ --physics.rangecut 1000000 \ """ @@ -29,6 +30,7 @@ rule beamline_acceptance_sim: --enableG4GPS \ --macroFile {input.macro} \ --compactFile $DETECTOR_PATH/epic_ip6_extended.xml \ + --printLevel WARNING \ --outputFile {output} \ --physics.rangecut 1000000 \ """ @@ -81,7 +83,7 @@ rule beamline_steering_reconstruction_preperation: script=workflow.source_path("processData.C"), data=SIMOUTDIR+"acceptanceTest{CAMPAIGN}.edm4hep.root", output: - onnxfile=ANALYSISDIR+"BeamlineSteeringReconstructionData{CAMPAIGN}.root", + rootfile=ANALYSISDIR+"BeamlineSteeringReconstructionData{CAMPAIGN}.root", shell: """ root -l -b -q '{input.script}("{input.data}", "{output.rootfile}")' @@ -102,13 +104,13 @@ rule beamline_steering_reconstruction_training: """ # Trains a regression model to predict the MCParticle momentum from the B2eR beamline virtual tracker hits. -rule beamline_steering_reconstruction_training: +rule beamline_steering_reconstruction_test: input: script=workflow.source_path("TestModel.py"), data=ANALYSISDIR+"BeamlineSteeringReconstructionData{CAMPAIGN}.root", model=ANALYSISDIR+"BeamlineSteeringReconstruction{CAMPAIGN}.onnx", output: - directory(ANALYSISDIR+"/NN_Test_{CAMPAIGN}/") + directory("sim_output/beamline/analysis/NN_Test_{CAMPAIGN}/") shell: """ python {input.script} --dataFiles {input.data} --modelFile {output.onnxfile} --outDir {output} diff --git a/benchmarks/beamline/acceptanceGPS.mac b/benchmarks/beamline/acceptanceGPS.mac index 59b48263..5f87735e 100644 --- a/benchmarks/beamline/acceptanceGPS.mac +++ b/benchmarks/beamline/acceptanceGPS.mac @@ -11,4 +11,4 @@ /gps/ene/min 6 GeV /gps/ene/max 18 GeV -/run/beamOn 1000000 \ No newline at end of file +/run/beamOn 4000000 \ No newline at end of file From 13030cca3e118ddd5510a72fddec1eea1c59d64d Mon Sep 17 00:00:00 2001 From: simonge Date: Wed, 25 Jun 2025 18:54:47 +0100 Subject: [PATCH 35/36] Fix snakemake attribute --- benchmarks/beamline/Snakefile | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/benchmarks/beamline/Snakefile b/benchmarks/beamline/Snakefile index 32338c9a..58294a94 100644 --- a/benchmarks/beamline/Snakefile +++ b/benchmarks/beamline/Snakefile @@ -113,7 +113,8 @@ rule beamline_steering_reconstruction_test: directory("sim_output/beamline/analysis/NN_Test_{CAMPAIGN}/") shell: """ - python {input.script} --dataFiles {input.data} --modelFile {output.onnxfile} --outDir {output} + mkdir {output} + python {input.script} --dataFiles {input.data} --modelFile {input.model} --outDir {output} """ rule beamline: From a208cd7d28312a21642574258c93aebaa0370c67 Mon Sep 17 00:00:00 2001 From: Simon Gardner Date: Mon, 7 Jul 2025 21:48:56 +0100 Subject: [PATCH 36/36] Update benchmarks/beamline/Snakefile Co-authored-by: Dmitry Kalinkin --- benchmarks/beamline/Snakefile | 1 + 1 file changed, 1 insertion(+) diff --git a/benchmarks/beamline/Snakefile b/benchmarks/beamline/Snakefile index 58294a94..75c524e7 100644 --- a/benchmarks/beamline/Snakefile +++ b/benchmarks/beamline/Snakefile @@ -20,6 +20,7 @@ rule beamline_steering_sim: rule beamline_acceptance_sim: input: + warmup="warmup/epic_ip6_extended.edm4hep.root", macro=workflow.source_path("acceptanceGPS.mac"), output: SIMOUTDIR+"acceptanceTest{CAMPAIGN}.edm4hep.root",