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

Update analysis

parent 841399b0
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
benchmarks/beamline/analysis.C
+ 341
− 32
View file @ 9a8c4d6d
...@@ -15,6 +15,7 @@ ...@@ -15,6 +15,7 @@
#include "TCanvas.h" #include "TCanvas.h"
#include "TStyle.h" #include "TStyle.h"
using RVecS = ROOT::VecOps::RVec<string>; using RVecS = ROOT::VecOps::RVec<string>;
using RNode = ROOT::RDF::RNode; using RNode = ROOT::RDF::RNode;
...@@ -23,12 +24,17 @@ void analysis( TString inFile = "/scratch/EIC/G4out/beamline/beamlineTest. ...@@ -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"){ std::string compactName = "/home/simong/EIC/epic/install/share/epic/epic_ip6_extended.xml"){
//Set ROOT style //Set ROOT style
gStyle->SetPadLeftMargin(0.25); // Set left margin gStyle->SetPadLeftMargin(0.1); // Set left margin
gStyle->SetPadRightMargin(0.15); // Set right margin gStyle->SetPadRightMargin(0.0); // Set right margin
gStyle->SetPadTopMargin(0.05); // Set top margin gStyle->SetPadTopMargin(0.0); // Set top margin
gStyle->SetPadBottomMargin(0.15);// Set bottom margin gStyle->SetPadBottomMargin(0.1);// Set bottom margin
gStyle->SetTitleYOffset(2); // Adjust y-axis title offset gStyle->SetTitleAlign(13);
gStyle->SetOptStat(1110); 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 //Set implicit multi-threading
ROOT::EnableImplicitMT(); ROOT::EnableImplicitMT();
...@@ -52,7 +58,37 @@ void analysis( TString inFile = "/scratch/EIC/G4out/beamline/beamlineTest. ...@@ -52,7 +58,37 @@ void analysis( TString inFile = "/scratch/EIC/G4out/beamline/beamlineTest.
if(Any(colNames==readoutName)){ if(Any(colNames==readoutName)){
d1 = d1.Define("pipeID",getSubID("pipe",detector),{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<volParams>& params) {
ROOT::VecOps::RVec<double> radii;
for (const auto& param : params) {
radii.push_back(param.radius);
}
return radii;
}, {"pipeParameters"})
.Define("xdet",[](const ROOT::VecOps::RVec<volParams>& params) {
ROOT::VecOps::RVec<double> xPos;
for (const auto& param : params) {
xPos.push_back(param.xPos);
}
return xPos;
}, {"pipeParameters"})
.Define("zdet",[](const ROOT::VecOps::RVec<volParams>& params) {
ROOT::VecOps::RVec<double> zPos;
for (const auto& param : params) {
zPos.push_back(param.zPos);
}
return zPos;
}, {"pipeParameters"})
.Define("rotation",[](const ROOT::VecOps::RVec<volParams>& params) {
ROOT::VecOps::RVec<double> rotation;
for (const auto& param : params) {
rotation.push_back(param.rotation);
}
return rotation;
}, {"pipeParameters"});
//global x,y,z position and momentum //global x,y,z position and momentum
d1 = d1.Define("xpos_global","BackwardsBeamlineHits.position.x") d1 = d1.Define("xpos_global","BackwardsBeamlineHits.position.x")
...@@ -66,9 +102,13 @@ void analysis( TString inFile = "/scratch/EIC/G4out/beamline/beamlineTest. ...@@ -66,9 +102,13 @@ void analysis( TString inFile = "/scratch/EIC/G4out/beamline/beamlineTest.
.Define("xpos","hitPosMom[0]") .Define("xpos","hitPosMom[0]")
.Define("ypos","hitPosMom[1]") .Define("ypos","hitPosMom[1]")
.Define("zpos","hitPosMom[2]") .Define("zpos","hitPosMom[2]")
.Define("xmom","hitPosMom[3]") .Define("xmomMag","hitPosMom[3]")
.Define("ymom","hitPosMom[4]") .Define("ymomMag","hitPosMom[4]")
.Define("zmom","hitPosMom[5]"); .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{ else{
...@@ -76,12 +116,57 @@ void analysis( TString inFile = "/scratch/EIC/G4out/beamline/beamlineTest. ...@@ -76,12 +116,57 @@ void analysis( TString inFile = "/scratch/EIC/G4out/beamline/beamlineTest.
return; 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; std::cout << "Executing Analysis and creating histograms" << std::endl;
//Create array of histogram results //Create array of histogram results
std::map<TString,ROOT::RDF::RResultPtr<TH2D>> hHistsxy;
std::map<TString,ROOT::RDF::RResultPtr<TH2D>> hHistsxyZoom;
std::map<TString,ROOT::RDF::RResultPtr<TH2D>> hHistsxpx; std::map<TString,ROOT::RDF::RResultPtr<TH2D>> hHistsxpx;
std::map<TString,ROOT::RDF::RResultPtr<TH2D>> hHistsypy; std::map<TString,ROOT::RDF::RResultPtr<TH2D>> hHistsypy;
std::map<TString,double> xMeans;
std::map<TString,double> yMeans;
std::map<TString,double> xStdDevs;
std::map<TString,double> yStdDevs;
std::map<TString,double> pxMeans;
std::map<TString,double> pyMeans;
std::map<TString,double> pxStdDevs;
std::map<TString,double> pyStdDevs;
//Fit paremeter and error maps
std::map<TString,double> xMeanFit;
std::map<TString,double> yMeanFit;
std::map<TString,double> xMeanFitErr;
std::map<TString,double> yMeanFitErr;
std::map<TString,double> xStdDevFit;
std::map<TString,double> yStdDevFit;
std::map<TString,double> xStdDevFitErr;
std::map<TString,double> yStdDevFitErr;
std::map<TString,double> pxMeanFit;
std::map<TString,double> pyMeanFit;
std::map<TString,double> pxMeanFitErr;
std::map<TString,double> pyMeanFitErr;
std::map<TString,double> pxStdDevFit;
std::map<TString,double> pyStdDevFit;
std::map<TString,double> pxStdDevFitErr;
std::map<TString,double> pyStdDevFitErr;
std::map<TString,double> pipeRadii;
std::map<TString,double> pipeXPos;
std::map<TString,double> pipeZPos;
std::map<TString,double> 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 //Create histograms
for(int i=0; i<=7; i++){ for(int i=0; i<=7; i++){
...@@ -91,54 +176,278 @@ void analysis( TString inFile = "/scratch/EIC/G4out/beamline/beamlineTest. ...@@ -91,54 +176,278 @@ void analysis( TString inFile = "/scratch/EIC/G4out/beamline/beamlineTest.
auto filterDF = d1.Define("xposf","xpos["+std::to_string(i)+"]") auto filterDF = d1.Define("xposf","xpos["+std::to_string(i)+"]")
.Define("yposf","ypos["+std::to_string(i)+"]") .Define("yposf","ypos["+std::to_string(i)+"]")
.Define("xmomf","xmom["+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 //Calculate Min and Max values
auto xmin = filterDF.Min("xposf").GetValue(); auto xminf = filterDF.Min("xposf").GetValue();
auto xmax = filterDF.Max("xposf").GetValue(); auto xmaxf = filterDF.Max("xposf").GetValue();
auto ymin = filterDF.Min("yposf").GetValue(); auto yminf = filterDF.Min("yposf").GetValue();
auto ymax = filterDF.Max("yposf").GetValue(); auto ymaxf = filterDF.Max("yposf").GetValue();
auto pxmin = filterDF.Min("xmomf").GetValue(); auto pxminf = filterDF.Min("xmomf").GetValue();
auto pxmax = filterDF.Max("xmomf").GetValue(); auto pxmaxf = filterDF.Max("xmomf").GetValue();
auto pymin = filterDF.Min("ymomf").GetValue(); auto pyminf = filterDF.Min("ymomf").GetValue();
auto pymax = filterDF.Max("ymomf").GetValue(); auto pymaxf = filterDF.Max("ymomf").GetValue();
// Calculate means and standard deviations
TString xname = name+";x offset [mm]; x trajectory component"; xMeans[name] = filterDF.Mean("xposf").GetValue();
TString yname = name+";y offset [mm]; y trajectory component"; yMeans[name] = filterDF.Mean("yposf").GetValue();
hHistsxpx[name] = filterDF.Histo2D({xname,xname,100,xmin,xmax,100,pxmin,pxmax},"xposf","xmomf"); xStdDevs[name] = filterDF.StdDev("xposf").GetValue();
hHistsypy[name] = filterDF.Histo2D({yname,yname,100,ymin,ymax,100,pymin,pymax},"yposf","ymomf"); 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); TCanvas *cX = new TCanvas("x_px_canvas","x_px_canvas",3000,1600);
cX->Divide(4,2); cX->Divide(4,2);
int i=1; i=1;
//Write histograms to file //Write histograms to file
for(auto& h : hHistsxpx){ for(auto& h : hHistsxpx){
cX->cd(i++); 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); TCanvas *cY = new TCanvas("y_py_canvas","y_py_canvas",3000,1600);
cY->Divide(4,2); cY->Divide(4,2);
i=1; i=1;
for(auto& h : hHistsypy){ for(auto& h : hHistsypy){
cY->cd(i++); 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"); TFile *f = new TFile(outFile,"RECREATE");
cXY->Write();
cX->Write(); cX->Write();
cY->Write(); cY->Write();
cFittedMeans->Write();
cFittedMomentumMeans->Write();
cPipeParams->Write();
f->Close(); f->Close();
std::cout << "Saving events to file" << std::endl; std::cout << "Saving events to file" << std::endl;
// ROOT::RDF::RSnapshotOptions opts; ROOT::RDF::RSnapshotOptions opts;
// opts.fMode = "UPDATE"; // 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); // 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);
} }
benchmarks/beamline/beamGPS.mac
+ 2
− 1
View file @ 9a8c4d6d
...@@ -4,7 +4,8 @@ ...@@ -4,7 +4,8 @@
/gps/pos/type Beam /gps/pos/type Beam
/gps/pos/sigma_x 0.119 mm /gps/pos/sigma_x 0.119 mm
/gps/pos/sigma_y 0.0107 mm /gps/pos/sigma_y 0.0107 mm
/gps/energy 18 GeV /gps/energy 17.846 GeV
#/gps/energy 18 GeV
/gps/particle e- /gps/particle e-
/run/beamOn 100000 /run/beamOn 100000
\ No newline at end of file
benchmarks/beamline/functors.h
+ 76
− 1
View file @ 9a8c4d6d
...@@ -101,3 +101,78 @@ struct globalToLocal{ ...@@ -101,3 +101,78 @@ struct globalToLocal{
dd4hep::Detector& detector; dd4hep::Detector& detector;
dd4hep::VolumeManager volumeManager; dd4hep::VolumeManager volumeManager;
}; };
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<volParams> operator()(const RVecHits& evt) {
ROOT::VecOps::RVec<volParams> 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<TString, double>& valueMap,
const std::map<TString, double>& 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
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment