diff --git a/benchmarks/beamline/Snakefile b/benchmarks/beamline/Snakefile
index 70eaf753da4eacc3fc841c92c3b110404dba3408..94b5c1e779fd1396b07d76eb05f1375f3b49e43d 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 a0551139c667e2dc9d069dc300e0f21d5c755adb..0c4c9465e9434d5e7f32974c928caccd35a6bbf7 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<string>;
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<TString,ROOT::RDF::RResultPtr<TH2D>> hHistsxpx;
+ std::map<TString,ROOT::RDF::RResultPtr<TH2D>> 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 d81b456765f689e8be2a0c63117c1dbe3d2ca79a..ee2414d08aa39c58919f7c4f795337a2c1f081bb 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<edm4hep::SimTrackerHitData>;
@@ -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<ROOT::VecOps::RVec<double>> operator()(const RVecHits& evt) {
+
+ ROOT::VecOps::RVec<ROOT::VecOps::RVec<double>> result;
+ ROOT::VecOps::RVec<double> xvec;
+ ROOT::VecOps::RVec<double> yvec;
+ ROOT::VecOps::RVec<double> zvec;
+ ROOT::VecOps::RVec<double> pxvec;
+ ROOT::VecOps::RVec<double> pyvec;
+ ROOT::VecOps::RVec<double> 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