benchmarks/beamline/beamlineAnalysis.C

0001 // Script to plot the x and y positions and momentum of beamline particles as they pass through the magnets
0002
0003 #include <iostream>
0004 #include "DD4hep/Detector.h"
0005 #include "DDRec/CellIDPositionConverter.h"
0006 #include "edm4hep/MCParticleCollection.h"
0007 #include "edm4hep/SimTrackerHitCollection.h"
0008 #include "edm4hep/SimCalorimeterHitCollection.h"
0009 #include "edm4hep/Vector3f.h"
0010 #include "edm4hep/Vector3d.h"
0011 #include "ROOT/RDataFrame.hxx"
0012 #include "ROOT/RDF/RInterface.hxx"
0013 #include "ROOT/RVec.hxx"
0014 #include "shared_functions.h"
0015 #include "TCanvas.h"
0016 #include "TStyle.h"
0017
0018
0019 using RVecS       = ROOT::VecOps::RVec<string>;
0020 using RNode       = ROOT::RDF::RNode;
0021
0022 int beamlineAnalysis(   TString inFile          = "/scratch/EIC/G4out/beamline/beamlineTest.edm4hep.root",
0023                         TString outFile         = "output.root",
0024                         std::string compactName = "/home/simong/EIC/epic/install/share/epic/epic_ip6_extended.xml",
0025                         TString beamspotCanvasName = "beamspot_canvas.png",
0026                         TString x_pxCanvasName = "x_px_canvas.png",
0027                         TString y_pyCanvasName = "y_py_canvas.png",
0028                         TString fittedPositionMeansCanvasName = "fitted_means_stddevs.png",
0029                         TString fittedMomentumMeansCanvasName = "fitted_momentum_means_stddevs.png",
0030                         TString pipeParamsCanvasName = "pipe_parameters.png"
0031                     ){
0032
0033     //Set ROOT style
0034     gStyle->SetPadLeftMargin(0.1);  // Set left margin
0035     gStyle->SetPadRightMargin(0.0); // Set right margin
0036     gStyle->SetPadTopMargin(0.0);   // Set top margin
0037     gStyle->SetPadBottomMargin(0.1);// Set bottom margin
0038     gStyle->SetTitleAlign(13);
0039     gStyle->SetTitleX(0.12);          // Place the title on the top right
0040     gStyle->SetTitleY(0.985);         // Place the title on the top right
0041     gStyle->SetTitleSize(0.08, "t");
0042     gStyle->SetTitleXOffset(1.0);    // Adjust y-axis title offset
0043     gStyle->SetTitleYOffset(1.0);    // Adjust y-axis title offset
0044     gStyle->SetOptStat(0);
0045
0046     //Set implicit multi-threading
0047     ROOT::EnableImplicitMT();
0048
0049     int pass = 0;
0050
0051     //Load the detector config
0052     dd4hep::Detector& detector = dd4hep::Detector::getInstance();
0053     detector.fromCompact(compactName);
0054
0055     ROOT::RDataFrame d0("events",inFile, {"BackwardsBeamlineHits"});
0056     RNode d1 = d0;
0057     RVecS colNames = d1.GetColumnNames();
0058
0059     //Set number of entries to process
0060     // d1 = d1.Range(0,1000);
0061     int   nEntries          = d1.Count().GetValue();
0062     float acceptableLoss    = 0.999; // Set the acceptable loss percentage to 0.1%
0063     float acceptableEntries = nEntries * acceptableLoss;
0064
0065     //Get the collection
0066     std::string readoutName = "BackwardsBeamlineHits";
0067
0068     std::cout << "Running lazy RDataframe execution" << std::endl;
0069
0070     if(Any(colNames==readoutName)){
0071
0072         d1 = d1.Define("pipeID",getSubID("pipe",detector),{readoutName})
0073                .Define("endID",getSubID("end",detector),{readoutName})
0074                .Define("pipeParameters",getVolumeParametersFromCellID(detector),{readoutName})
0075                .Define("pipeRadius",[](const ROOT::VecOps::RVec<volParams>& params) {
0076                 ROOT::VecOps::RVec<double> radii;
0077                 for (const auto& param : params) {
0078                     radii.push_back(param.radius);
0079                 }
0080                 return radii;
0081                 }, {"pipeParameters"})
0082                 .Define("xdet",[](const ROOT::VecOps::RVec<volParams>& params) {
0083                 ROOT::VecOps::RVec<double> xPos;
0084                 for (const auto& param : params) {
0085                     xPos.push_back(param.xPos);
0086                 }
0087                 return xPos;
0088                 }, {"pipeParameters"})
0089                 .Define("zdet",[](const ROOT::VecOps::RVec<volParams>& params) {
0090                 ROOT::VecOps::RVec<double> zPos;
0091                 for (const auto& param : params) {
0092                     zPos.push_back(param.zPos);
0093                 }
0094                 return zPos;
0095                 }, {"pipeParameters"})
0096                 .Define("rotation",[](const ROOT::VecOps::RVec<volParams>& params) {
0097                 ROOT::VecOps::RVec<double> rotation;
0098                 for (const auto& param : params) {
0099                     rotation.push_back(param.rotation);
0100                 }
0101                 return rotation;
0102                 }, {"pipeParameters"});
0103
0104
0105         //global x,y,z position and momentum
0106         d1 = d1.Define("xpos_global","BackwardsBeamlineHits.position.x")
0107                 .Define("ypos_global","BackwardsBeamlineHits.position.y")
0108                 .Define("zpos_global","BackwardsBeamlineHits.position.z")
0109                 .Define("px_global","BackwardsBeamlineHits.momentum.x")
0110                 .Define("py_global","BackwardsBeamlineHits.momentum.y")
0111                 .Define("pz_global","BackwardsBeamlineHits.momentum.z");
0112
0113         d1 = d1.Define("hitPosMom",globalToLocal(detector),{readoutName})
0114                 .Define("xpos","hitPosMom[0]")
0115                 .Define("ypos","hitPosMom[1]")
0116                 .Define("zpos","hitPosMom[2]")
0117                 .Define("xmomMag","hitPosMom[3]")
0118                 .Define("ymomMag","hitPosMom[4]")
0119                 .Define("zmomMag","hitPosMom[5]")
0120                 .Define("momMag","sqrt(xmomMag*xmomMag+ymomMag*ymomMag+zmomMag*zmomMag)")
0121                 .Define("xmom","xmomMag/momMag")
0122                 .Define("ymom","ymomMag/momMag")
0123                 .Define("zmom","zmomMag/momMag");
0124
0125     }
0126     else{
0127         std::cout << "Collection " << readoutName << " not found in file" << std::endl;
0128         return 1;
0129     }
0130
0131     // Calculate the maximum pipe radius for plotting
0132     auto maxPipeRadius = 1.2*d1.Max("pipeRadius").GetValue();
0133
0134     std::cout << "Executing Analysis and creating histograms" << std::endl;
0135
0136     //Create array of histogram results
0137     std::map<TString,ROOT::RDF::RResultPtr<TH2D>> hHistsxy;
0138     std::map<TString,ROOT::RDF::RResultPtr<TH2D>> hHistsxyZoom;
0139     std::map<TString,ROOT::RDF::RResultPtr<TH2D>> hHistsxpx;
0140     std::map<TString,ROOT::RDF::RResultPtr<TH2D>> hHistsypy;
0141
0142     std::map<TString,double> xMeans;
0143     std::map<TString,double> yMeans;
0144     std::map<TString,double> xStdDevs;
0145     std::map<TString,double> yStdDevs;
0146     std::map<TString,double> pxMeans;
0147     std::map<TString,double> pyMeans;
0148     std::map<TString,double> pxStdDevs;
0149     std::map<TString,double> pyStdDevs;
0150
0151     //Fit paremeter and error maps
0152     std::map<TString,double> xMeanFit;
0153     std::map<TString,double> yMeanFit;
0154     std::map<TString,double> xMeanFitErr;
0155     std::map<TString,double> yMeanFitErr;
0156     std::map<TString,double> xStdDevFit;
0157     std::map<TString,double> yStdDevFit;
0158     std::map<TString,double> xStdDevFitErr;
0159     std::map<TString,double> yStdDevFitErr;
0160     std::map<TString,double> pxMeanFit;
0161     std::map<TString,double> pyMeanFit;
0162     std::map<TString,double> pxMeanFitErr;
0163     std::map<TString,double> pyMeanFitErr;
0164     std::map<TString,double> pxStdDevFit;
0165     std::map<TString,double> pyStdDevFit;
0166     std::map<TString,double> pxStdDevFitErr;
0167     std::map<TString,double> pyStdDevFitErr;
0168
0169     std::map<TString,double> pipeRadii;
0170     std::map<TString,double> pipeXPos;
0171     std::map<TString,double> pipeZPos;
0172     std::map<TString,double> pipeRotation;
0173
0174     // Queue up all actions
0175     auto xmin_ptr  = d1.Min("xpos");
0176     auto xmax_ptr  = d1.Max("xpos");
0177     auto ymin_ptr  = d1.Min("ypos");
0178     auto ymax_ptr  = d1.Max("ypos");
0179     auto pxmin_ptr = d1.Min("xmom");
0180     auto pxmax_ptr = d1.Max("xmom");
0181     auto pymin_ptr = d1.Min("ymom");
0182     auto pymax_ptr = d1.Max("ymom");
0183
0184     // Now trigger the event loop (only once)
0185     auto xmin  = xmin_ptr.GetValue();
0186     auto xmax  = xmax_ptr.GetValue();
0187     auto ymin  = ymin_ptr.GetValue();
0188     auto ymax  = ymax_ptr.GetValue();
0189     auto pxmin = pxmin_ptr.GetValue();
0190     auto pxmax = pxmax_ptr.GetValue();
0191     auto pymin = pymin_ptr.GetValue();
0192     auto pymax = pymax_ptr.GetValue();
0193
0194     //Create histograms
0195     for(int i=0; i<=7; i++){
0196
0197         std::string name = "pipeID";
0198         name += std::to_string(i);
0199         auto filterDF = d1.Define("xposf","xpos[pipeID=="+std::to_string(i)+"]")
0200                           .Define("yposf","ypos[pipeID=="+std::to_string(i)+"]")
0201                           .Define("xmomf","xmom[pipeID=="+std::to_string(i)+"]")
0202                           .Define("ymomf","ymom[pipeID=="+std::to_string(i)+"]")
0203                           .Define("pipeRadiusf","pipeRadius[pipeID=="+std::to_string(i)+"]")
0204                           .Define("xdetf","xdet[pipeID=="+std::to_string(i)+"]")
0205                           .Define("zdetf","zdet[pipeID=="+std::to_string(i)+"]")
0206                           .Define("rotationf","rotation[pipeID=="+std::to_string(i)+"]");
0207
0208
0209         //Calculate Min and Max values
0210         auto xminf = filterDF.Min("xposf").GetValue();
0211         auto xmaxf = filterDF.Max("xposf").GetValue();
0212         auto yminf = filterDF.Min("yposf").GetValue();
0213         auto ymaxf = filterDF.Max("yposf").GetValue();
0214         auto pxminf = filterDF.Min("xmomf").GetValue();
0215         auto pxmaxf = filterDF.Max("xmomf").GetValue();
0216         auto pyminf = filterDF.Min("ymomf").GetValue();
0217         auto pymaxf = filterDF.Max("ymomf").GetValue();
0218         // If the min and max values are equal, set them to a small value
0219         if(xminf == xmaxf) {
0220             xminf -= 0.01;
0221             xmaxf += 0.01;
0222             pass = 1; // Set pass to 1 if xminf == xmaxf
0223             std::cout << "Warning: xminf == xmaxf for pipe ID " << i << ", likely no hits." << std::endl;
0224         }
0225         if(yminf == ymaxf) {
0226             yminf -= 0.01;
0227             ymaxf += 0.01;
0228             pass = 1; // Set pass to 1 if yminf == ymaxf
0229             std::cout << "Warning: yminf == ymaxf for pipe ID " << i << ", likely no hits." << std::endl;
0230         }
0231         if(pxminf == pxmaxf) {
0232             pxminf -= 0.01;
0233             pxmaxf += 0.01;
0234             pass = 1; // Set pass to 1 if pxminf == pxmaxf
0235             std::cout << "Warning: pxminf == pxmaxf for pipe ID " << i << ", likely no hits." << std::endl;
0236         }
0237         if(pyminf == pymaxf) {
0238             pyminf -= 0.01;
0239             pymaxf += 0.01;
0240             pass = 1; // Set pass to 1 if pyminf == pymaxf
0241             std::cout << "Warning: pyminf == pymaxf for pipe ID " << i << ", likely no hits." << std::endl;
0242         }
0243
0244         // Calculate means and standard deviations
0245         xMeans[name] = filterDF.Mean("xposf").GetValue();
0246         yMeans[name] = filterDF.Mean("yposf").GetValue();
0247         xStdDevs[name] = filterDF.StdDev("xposf").GetValue();
0248         yStdDevs[name] = filterDF.StdDev("yposf").GetValue();
0249         pxMeans[name] = filterDF.Mean("xmomf").GetValue();
0250         pyMeans[name] = filterDF.Mean("ymomf").GetValue();
0251         pxStdDevs[name] = filterDF.StdDev("xmomf").GetValue();
0252         pyStdDevs[name] = filterDF.StdDev("ymomf").GetValue();
0253
0254         // Calculate axes for zoomed beamspot histogram so that it is quare around the mean x and y
0255         double halfrange = std::max({xMeans[name]-xminf, xmaxf-xMeans[name], yMeans[name]-yminf, ymaxf-yMeans[name]});
0256         double xMinZoom = xMeans[name] - halfrange;
0257         double xMaxZoom = xMeans[name] + halfrange;
0258         double yMinZoom = yMeans[name] - halfrange;
0259         double yMaxZoom = yMeans[name] + halfrange;
0260
0261         TString beamspotName = "Beamspot ID"+std::to_string(i)+";x offset [cm]; y offset [cm]";
0262         TString xyname = name+";x Offset [cm]; y Offset [cm]";
0263         TString xname = name+";x Offset [cm]; x trajectory component";
0264         TString yname = name+";y Offset [cm]; y trajectory component";
0265         hHistsxy[name] = filterDF.Histo2D({beamspotName,beamspotName,400,-maxPipeRadius,maxPipeRadius,400,-maxPipeRadius,maxPipeRadius},"xposf","yposf");
0266         hHistsxyZoom[name] = filterDF.Histo2D({xyname,xyname,100,xMinZoom,xMaxZoom,100,yMinZoom,yMaxZoom},"xposf","yposf");
0267         hHistsxpx[name]    = filterDF.Histo2D({xname,xname,400,xmin,xmax,400,pxmin,pxmax},"xposf","xmomf");
0268         hHistsypy[name]    = filterDF.Histo2D({yname,yname,400,ymin,ymax,400,pymin,pymax},"yposf","ymomf");
0269
0270         //Parameters of the pipe
0271         pipeRadii[name]    = filterDF.Max("pipeRadiusf").GetValue();
0272         std::cout << "Pipe ID: " << name << " Radius: " << pipeRadii[name] << std::endl;
0273         pipeXPos[name]     = filterDF.Max("xdetf").GetValue();
0274         pipeZPos[name]     = filterDF.Max("zdetf").GetValue();
0275         pipeRotation[name] = filterDF.Max("rotationf").GetValue();
0276
0277         //Fit gaussian to the x, y, px and py distributions
0278         auto xhist = hHistsxy[name]->ProjectionX();
0279         auto yhist = hHistsxy[name]->ProjectionY();
0280         auto pxhist = hHistsxpx[name]->ProjectionY();
0281         auto pyhist = hHistsypy[name]->ProjectionY();
0282         xhist->Fit("gaus","Q");
0283         yhist->Fit("gaus","Q");
0284         pxhist->Fit("gaus","Q");
0285         pyhist->Fit("gaus","Q");
0286         //Get the fit parameters and errors
0287         xMeanFit[name] = xhist->GetFunction("gaus")->GetParameter(1);
0288         yMeanFit[name] = yhist->GetFunction("gaus")->GetParameter(1);
0289         xMeanFitErr[name] = xhist->GetFunction("gaus")->GetParError(1);
0290         yMeanFitErr[name] = yhist->GetFunction("gaus")->GetParError(1);
0291         xStdDevFit[name] = xhist->GetFunction("gaus")->GetParameter(2);
0292         yStdDevFit[name] = yhist->GetFunction("gaus")->GetParameter(2);
0293         xStdDevFitErr[name] = xhist->GetFunction("gaus")->GetParError(2);
0294         yStdDevFitErr[name] = yhist->GetFunction("gaus")->GetParError(2);
0295         pxMeanFit[name] = pxhist->GetFunction("gaus")->GetParameter(1);
0296         pyMeanFit[name] = pyhist->GetFunction("gaus")->GetParameter(1);
0297         pxMeanFitErr[name] = pxhist->GetFunction("gaus")->GetParError(1);
0298         pyMeanFitErr[name] = pyhist->GetFunction("gaus")->GetParError(1);
0299         pxStdDevFit[name] = pxhist->GetFunction("gaus")->GetParameter(2);
0300         pyStdDevFit[name] = pyhist->GetFunction("gaus")->GetParameter(2);
0301         pxStdDevFitErr[name] = pxhist->GetFunction("gaus")->GetParError(2);
0302         pyStdDevFitErr[name] = pyhist->GetFunction("gaus")->GetParError(2);
0303
0304     }
0305
0306
0307
0308     // Create histograms of the beamspot
0309     TCanvas *cXY = new TCanvas("beamspot_canvas","beamspot_canvas",3000,1600);
0310     cXY->Divide(4,2);
0311     int i=1;
0312     for(auto [name,h] : hHistsxy){
0313
0314         // Check histogram entries
0315         if(h->GetEntries() < acceptableEntries){
0316             std::cout << "Warning: Only " << h->GetEntries()/nEntries << " of particles contributing to histogram " <<  name
0317                       << " , which is below the accepted threshold of " << acceptableEntries/nEntries << std::endl;
0318             pass = 1;
0319         }
0320
0321         // Get the pipe radius for this histogram
0322         auto pipeRadius = pipeRadii[name];
0323         cXY->cd(i++);
0324
0325         h->Draw("col");
0326         //Draw cicle
0327         TEllipse *circle = new TEllipse(0,0,pipeRadius);
0328         circle->SetLineColor(kRed);
0329         circle->SetLineWidth(2);
0330         circle->SetFillStyle(0);
0331         circle->Draw("same");
0332
0333         // Add zoomed version in the top-right corner
0334         TPad *pad = new TPad("zoomPad", "Zoomed View", 0.65, 0.65, 1.0, 1.0);
0335         pad->SetFillStyle(0); // Transparent background
0336         pad->Draw();
0337         pad->cd();
0338         // Draw the zoomed histogram without its title or axis titles
0339         hHistsxyZoom[name]->SetTitle("");
0340         hHistsxyZoom[name]->GetXaxis()->SetTitle("");
0341         hHistsxyZoom[name]->GetYaxis()->SetTitle("");
0342         hHistsxyZoom[name]->Draw("col");
0343         cXY->cd(); // Return to the main canvas
0344     }
0345
0346     // x-px canvas
0347     TCanvas *cX = new TCanvas("x_px_canvas","x_px_canvas",3000,1600);
0348     cX->Divide(4,2);
0349
0350     i=1;
0351     //Write histograms to file
0352     for(auto& h : hHistsxpx){
0353         cX->cd(i++);
0354         h.second->Draw("col");
0355     }
0356
0357     // y-py canvas
0358     TCanvas *cY = new TCanvas("y_py_canvas","y_py_canvas",3000,1600);
0359     cY->Divide(4,2);
0360
0361     i=1;
0362     for(auto& h : hHistsypy){
0363         cY->cd(i++);
0364         h.second->Draw("col");
0365     }
0366
0367     // Save 2D canvases as pngs
0368     cXY->SaveAs(beamspotCanvasName);
0369     cX->SaveAs(x_pxCanvasName);
0370     cY->SaveAs(y_pyCanvasName);
0371
0372     // ---------------------------------------------------------------------------
0373     // Create histograms showing the fitted means and standard deviations of the positions and momenta
0374     // ---------------------------------------------------------------------------
0375
0376     // Create histograms for fitted X means and standard deviations
0377     TH1F* hFittedXMeans = CreateFittedHistogram("hFittedXMeans",
0378         "Mean X Offset [cm]",
0379         xMeanFit,
0380         xMeanFitErr,
0381         "Pipe ID");
0382
0383     TH1F* hFittedXStdDevs = CreateFittedHistogram("hFittedXStdDevs",
0384         "Std Deviation X Offset [cm]",
0385         xStdDevFit,
0386         xStdDevFitErr,
0387         "Pipe ID");
0388
0389     // Create histograms for fitted Y means and standard deviations
0390     TH1F* hFittedYMeans = CreateFittedHistogram("hFittedYMeans",
0391         "Mean Y Offset [cm]",
0392         yMeanFit,
0393         yMeanFitErr,
0394         "Pipe ID");
0395
0396     TH1F* hFittedYStdDevs = CreateFittedHistogram("hFittedYStdDevs",
0397         "Std Deviation Y Offset [cm]",
0398         yStdDevFit,
0399         yStdDevFitErr,
0400         "Pipe ID");
0401
0402     TH1F* hFittedPxMeans = CreateFittedHistogram("hFittedPxMeans",
0403         "Mean Px",
0404         pxMeanFit,
0405         pxMeanFitErr,
0406         "Pipe ID");
0407     TH1F* hFittedPyMeans = CreateFittedHistogram("hFittedPyMeans",
0408         "Mean Py",
0409         pyMeanFit,
0410         pyMeanFitErr,
0411         "Pipe ID");
0412     TH1F* hFittedPxStdDevs = CreateFittedHistogram("hFittedPxStdDevs",
0413         "Std Deviation Px",
0414         pxStdDevFit,
0415         pxStdDevFitErr,
0416         "Pipe ID");
0417     TH1F* hFittedPyStdDevs = CreateFittedHistogram("hFittedPyStdDevs",
0418         "Std Deviation Py",
0419         pyStdDevFit,
0420         pyStdDevFitErr,
0421         "Pipe ID");
0422
0423
0424     // Create a canvas for the fitted beamspot means and standard deviations
0425     TCanvas *cFittedMeans = new TCanvas("cFittedMeans", "Fitted Beamspot Means and Std Deviation", 1200, 800);
0426     cFittedMeans->Divide(2, 2);
0427     cFittedMeans->cd(1);
0428     hFittedXMeans->Draw("E1"); // "E1" draws error bars
0429     cFittedMeans->cd(2);
0430     hFittedXStdDevs->Draw("E1"); // "E1" draws error bars
0431     cFittedMeans->cd(3);
0432     hFittedYMeans->Draw("E1"); // "E1" draws error bars
0433     cFittedMeans->cd(4);
0434     hFittedYStdDevs->Draw("E1"); // "E1" draws error bars
0435     cFittedMeans->SetGrid();
0436     cFittedMeans->Update();
0437     // Save the canvas as a PNG file
0438     cFittedMeans->SaveAs(fittedPositionMeansCanvasName);
0439
0440     // Create a canvas for the fitted momentum means and standard deviations
0441     TCanvas *cFittedMomentumMeans = new TCanvas("cFittedMomentumMeans", "Fitted Momentum Means and Std Deviation", 1200, 800);
0442     cFittedMomentumMeans->Divide(2, 2);
0443     cFittedMomentumMeans->cd(1);
0444     hFittedPxMeans->Draw("E1"); // "E1" draws error bars
0445     cFittedMomentumMeans->cd(2);
0446     hFittedPxStdDevs->Draw("E1"); // "E1" draws error bars
0447     cFittedMomentumMeans->cd(3);
0448     hFittedPyMeans->Draw("E1"); // "E1" draws error bars
0449     cFittedMomentumMeans->cd(4);
0450     hFittedPyStdDevs->Draw("E1"); // "E1" draws error bars
0451     cFittedMomentumMeans->SetGrid();
0452     cFittedMomentumMeans->Update();
0453     // Save the canvas as a PNG file
0454     cFittedMomentumMeans->SaveAs(fittedMomentumMeansCanvasName);
0455
0456
0457     // -----------------------------------------------------------------------------
0458     // Create histograms of the beampipe parameters
0459     // -----------------------------------------------------------------------------
0460     TH1F* hPipeRadii = CreateFittedHistogram("hPipeRadii",
0461         "Radius [cm]",
0462         pipeRadii,
0463         {},
0464         "Pipe ID");
0465     TH1F* hPipeXPos = CreateFittedHistogram("hPipeXPos",
0466         "X Position [cm]",
0467         pipeXPos,
0468         {},
0469         "Pipe ID");
0470     TH1F* hPipeZPos = CreateFittedHistogram("hPipeZPos",
0471         "Z Position [cm]",
0472         pipeZPos,
0473         {},
0474         "Pipe ID");
0475     TH1F* hPipeRotations = CreateFittedHistogram("hPipeRotations",
0476         "Rotation [rad]",
0477         pipeRotation,
0478         {},
0479         "Pipe ID");
0480
0481     // Create a canvas for the pipe parameters
0482     TCanvas *cPipeParams = new TCanvas("cPipeParams", "Pipe Parameters", 1200, 400);
0483     cPipeParams->Divide(4, 1);
0484     cPipeParams->cd(1);
0485     hPipeRadii->Draw("");
0486     cPipeParams->cd(2);
0487     hPipeXPos->Draw("");
0488     cPipeParams->cd(3);
0489     hPipeZPos->Draw("");
0490     cPipeParams->cd(4);
0491     hPipeRotations->Draw("");
0492     cPipeParams->SetGrid();
0493     cPipeParams->Update();
0494     // Save the canvas as a PNG file
0495     cPipeParams->SaveAs(pipeParamsCanvasName);
0496
0497
0498     TFile *f = new TFile(outFile,"RECREATE");
0499     cXY->Write();
0500     cX->Write();
0501     cY->Write();
0502     cFittedMeans->Write();
0503     cFittedMomentumMeans->Write();
0504     cPipeParams->Write();
0505
0506     f->Close();
0507
0508     std::cout << "Analysis complete. Results saved to " << outFile << std::endl;
0509     return pass;
0510
0511     // std::cout << "Saving events to file" << std::endl;
0512
0513     // ROOT::RDF::RSnapshotOptions opts;
0514     //  opts.fMode = "UPDATE";
0515     //  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);
0516 }