Skip to content
Snippets Groups Projects
Commit 5b476975 authored by Simon Gardner's avatar Simon Gardner
Browse files

Update

parent a853f58e
No related branches found
No related tags found
No related merge requests found
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}")'
"""
......@@ -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);
}
#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
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment