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