Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-12-16 09:27:47 

0001 // Plots the resolution of the reconstructed particles
0002 
0003 #include <iostream>
0004 #include "ROOT/RDataFrame.hxx"
0005 #include "ROOT/RDF/RInterface.hxx"
0006 #include "ROOT/RVec.hxx"
0007 #include "edm4hep/MCParticleCollection.h"
0008 #include "edm4eic/ReconstructedParticleCollection.h"
0009 #include "TCanvas.h"
0010 #include "TFile.h"
0011 #include "TStyle.h"
0012 
0013 void reconstructionAnalysis(TString inFile             = "/home/simong/EIC/detector_benchmarks_anl/sim_output/lowq2_reconstruction/analysis/Low-Q2_retrained_Particles_new.eicrecon.edm4hep.root",
0014                             float   beamEnergy         = 18.0,
0015                             TString outFile            = "reconstruction_results.root",
0016                             TString momentumCanvasName = "momentum_resolution.png",
0017                             TString energyThetaPhiCanvasName = "energy_theta_phi_resolution.png",
0018                             TString relationCanvasName = "relation_resolution.png",
0019                             TString resolutionGraphsCanvasName = "resolution_graphs.png",
0020                             std::string particleCollectionName = "TaggerTrackerReconstructedParticles") {
0021 
0022     //Set ROOT style    
0023     gStyle->SetPadLeftMargin(0.1);  // Set left margin
0024     gStyle->SetPadRightMargin(0.0); // Set right margin
0025     gStyle->SetPadTopMargin(0.0);   // Set top margin
0026     gStyle->SetPadBottomMargin(0.1);// Set bottom margin
0027     gStyle->SetTitleAlign(13);
0028     gStyle->SetTitleX(0.12);          // Place the title on the top right
0029     gStyle->SetTitleY(0.985);         // Place the title on the top right
0030     gStyle->SetTitleSize(0.08, "t");
0031     gStyle->SetTitleXOffset(1.0);    // Adjust y-axis title offset
0032     gStyle->SetTitleYOffset(1.0);    // Adjust y-axis title offset
0033     gStyle->SetOptStat(0);
0034 
0035     ROOT::RDataFrame d0("events",inFile, {"MCParticles",particleCollectionName});
0036 
0037     auto filterDF = d0.Define("SimParticles", "MCParticles[MCParticles.generatorStatus==1 && MCParticles.PDG==11]")
0038                       .Filter("SimParticles.size()==1")
0039                       .Filter(particleCollectionName+".size()==1");
0040 
0041     // Plot x,y,z momentum resolution as a function of the MCParticle momentum
0042     auto momentumDF = filterDF
0043         .Define("reco_momentum", particleCollectionName+"[0].momentum")
0044         .Define("mc_momentum", "SimParticles[0].momentum")
0045         .Define("px_rec", "reco_momentum.x")
0046         .Define("py_rec", "reco_momentum.y")
0047         .Define("pz_rec", "reco_momentum.z")
0048         .Define("px_mc", "mc_momentum.x")
0049         .Define("py_mc", "mc_momentum.y")
0050         .Define("pz_mc", "mc_momentum.z")
0051         // Calculate theta and phi for reco and MC
0052         .Define("theta_rec", "std::atan2(std::sqrt(px_rec*px_rec + py_rec*py_rec), pz_rec)")
0053         .Define("theta_mc", "std::atan2(std::sqrt(px_mc*px_mc + py_mc*py_mc), pz_mc)")
0054         .Define("phi_rec_rad", "std::atan2(py_rec, px_rec)")
0055         .Define("phi_mc_rad", "std::atan2(py_mc, px_mc)")
0056         .Define("phi_rec", "TMath::RadToDeg()*phi_rec_rad")
0057         .Define("phi_mc", "TMath::RadToDeg()*phi_mc_rad")
0058         // Calculate resolutions
0059         .Define("px_diff", "(px_rec - px_mc)")
0060         .Define("py_diff", "(py_rec - py_mc)")
0061         .Define("pz_res", "(pz_rec - pz_mc)/pz_mc")
0062         .Define("E_rec", particleCollectionName+"[0].energy")
0063         .Define("E_mc", "std::sqrt(px_mc*px_mc + py_mc*py_mc + pz_mc*pz_mc + SimParticles[0].mass*SimParticles[0].mass)")
0064         .Define("E_res", "(E_rec - E_mc)/E_mc")
0065         .Define("theta_diff", "(theta_rec - theta_mc)")
0066         .Define("phi_diff", "TMath::RadToDeg()*ROOT::VecOps::DeltaPhi(phi_rec_rad, phi_mc_rad)");
0067 
0068     //Print the size of the original DataFrame
0069     std::cout << "Original DataFrame size: " << d0.Count().GetValue() << std::endl;
0070     //Print the size of the filtered DataFrame
0071     std::cout << "Filtered DataFrame size: " << filterDF.Count().GetValue() << std::endl;
0072 
0073     int   momentumBins      = 100;
0074     float momentumXRange[2] = {-0.1, 0.1};
0075     float momentumYRange[2] = {-0.1, 0.1};
0076     float momentumZRange[2] = {-beamEnergy, 0};
0077 
0078     int   momentumResolutionBins      = 100;
0079     float momentumDifferenceXRange[2] = {-0.05, 0.05};
0080     float momentumDifferenceYRange[2] = {-0.05, 0.05};
0081     float momentumResolutionZRange[2] = {-0.1, 0.1};
0082 
0083     int   energyBins      = 100;
0084     int   thetaBins       = 100;
0085     int   phiBins         = 100;
0086     float energyRange[2]  = {3.5, beamEnergy}; // GeV
0087     float thetaRange[2]   = {3.134, TMath::Pi()}; // radians from 3.1 to pi
0088     float phiRange[2]     = {-180, 180}; // degrees from -180 to 180
0089 
0090     int   resolutionBins           = 100;
0091     float energyResolutionRange[2] = {-0.1, 0.1};
0092     float thetaResolutionRange[2]  = {-0.003, 0.003};
0093     float phiResolutionRange[2]    = {-90, 90}; // degrees from -90 to 90
0094 
0095     // Plot reconstructed vs montecarlo momentum components
0096     auto px_Hist = momentumDF.Histo2D({"px_vs_px", "Reconstructed vs MC px; px reconstructed [GeV]; px MC [GeV]", momentumBins, momentumXRange[0], momentumXRange[1], momentumBins, momentumXRange[0], momentumXRange[1]}, "px_rec", "px_mc");
0097     auto py_Hist = momentumDF.Histo2D({"py_vs_py", "Reconstructed vs MC py; py reconstructed [GeV]; py MC [GeV]", momentumBins, momentumYRange[0], momentumYRange[1], momentumBins, momentumYRange[0], momentumYRange[1]}, "py_rec", "py_mc");
0098     auto pz_Hist = momentumDF.Histo2D({"pz_vs_pz", "Reconstructed vs MC pz; pz reconstructed [GeV]; pz MC [GeV]", momentumBins, momentumZRange[0], momentumZRange[1], momentumBins, momentumZRange[0], momentumZRange[1]}, "pz_rec", "pz_mc");
0099 
0100     // Plot individual momentum resolutions
0101     auto px_diff_Hist = momentumDF.Histo1D({"px_diff", "px difference; px difference [GeV]; Entries", momentumResolutionBins, momentumDifferenceXRange[0], momentumDifferenceXRange[1]}, "px_diff");
0102     auto py_diff_Hist = momentumDF.Histo1D({"py_diff", "py difference; py difference [GeV]; Entries", momentumResolutionBins, momentumDifferenceYRange[0], momentumDifferenceYRange[1]}, "py_diff");
0103     auto pz_res_Hist  = momentumDF.Histo1D({"pz_res",  "pz resolution; pz resolution [GeV]; Entries", momentumResolutionBins, momentumResolutionZRange[0], momentumResolutionZRange[1]}, "pz_res");
0104 
0105     // Plot reconstructed vs montecarlo energy, theta and phi
0106     auto E_Hist     = momentumDF.Histo2D({"E_vs_E",         "Reconstructed vs MC energy; E reconstructed [GeV]; E MC [GeV]",        energyBins, energyRange[0], energyRange[1], energyBins, energyRange[0], energyRange[1]}, "E_rec", "E_mc");
0107     auto theta_Hist = momentumDF.Histo2D({"theta_vs_theta", "Reconstructed vs MC theta; theta reconstructed [rad]; theta MC [rad]", thetaBins, thetaRange[0], thetaRange[1], thetaBins, thetaRange[0], thetaRange[1]}, "theta_rec", "theta_mc");
0108     auto phi_Hist   = momentumDF.Histo2D({"phi_vs_phi",     "Reconstructed vs MC phi; phi reconstructed [deg]; phi MC [deg]",       phiBins, phiRange[0], phiRange[1], phiBins, phiRange[0], phiRange[1]}, "phi_rec", "phi_mc");
0109 
0110     auto E_res_Hist      = momentumDF.Histo1D({"E_res", "E resolution; E resolution [GeV]; Entries", resolutionBins, energyResolutionRange[0], energyResolutionRange[1]}, "E_res");
0111     auto theta_diff_Hist = momentumDF.Histo1D({"theta_diff", "theta difference; theta difference [rad]; Entries", resolutionBins, thetaResolutionRange[0], thetaResolutionRange[1]}, "theta_diff");
0112     auto phi_diff_Hist   = momentumDF.Histo1D({"phi_diff", "phi difference; phi difference [deg]; Entries", resolutionBins, phiResolutionRange[0], phiResolutionRange[1]}, "phi_diff");
0113 
0114     // Plot Reconstructed energy, theta and phi resolutions as a function of each reconstructed value of energy, thata and phi
0115     auto E_res_vs_E_Hist          = momentumDF.Histo2D({"E_res_vs_E", "E resolution vs E reconstructed; E reconstructed [GeV]; E resolution [GeV]", energyBins, energyRange[0], energyRange[1], resolutionBins, energyResolutionRange[0], energyResolutionRange[1]}, "E_rec", "E_res");
0116     auto E_res_vs_theta_Hist      = momentumDF.Histo2D({"E_res_vs_theta", "E resolution vs theta reconstructed; theta reconstructed [rad]; E resolution [GeV]", thetaBins, thetaRange[0], thetaRange[1], resolutionBins, energyResolutionRange[0], energyResolutionRange[1]}, "theta_rec", "E_res");
0117     auto E_res_vs_phi_Hist        = momentumDF.Histo2D({"E_res_vs_phi", "E resolution vs phi reconstructed; phi reconstructed [deg]; E resolution [GeV]", phiBins, phiRange[0], phiRange[1], resolutionBins, energyResolutionRange[0], energyResolutionRange[1]}, "phi_rec", "E_res");
0118     auto theta_diff_vs_E_Hist     = momentumDF.Histo2D({"theta_diff_vs_E", "theta difference vs E reconstructed; E reconstructed [GeV]; theta difference [rad]", energyBins, energyRange[0], energyRange[1], resolutionBins, thetaResolutionRange[0], thetaResolutionRange[1]}, "E_rec", "theta_diff");
0119     auto theta_diff_vs_theta_Hist = momentumDF.Histo2D({"theta_diff_vs_theta", "theta difference vs theta reconstructed; theta reconstructed [rad]; theta difference [rad]", thetaBins, thetaRange[0], thetaRange[1], resolutionBins, thetaResolutionRange[0], thetaResolutionRange[1]}, "theta_rec", "theta_diff");
0120     auto theta_diff_vs_phi_Hist   = momentumDF.Histo2D({"theta_diff_vs_phi", "theta difference vs phi reconstructed; phi reconstructed [deg]; theta difference [rad]", phiBins, phiRange[0], phiRange[1], resolutionBins, thetaResolutionRange[0], thetaResolutionRange[1]}, "phi_rec", "theta_diff");
0121     auto phi_diff_vs_E_Hist       = momentumDF.Histo2D({"phi_diff_vs_E", "phi difference vs E reconstructed; E reconstructed [GeV]; phi difference [rad]", energyBins, energyRange[0], energyRange[1], resolutionBins, phiResolutionRange[0], phiResolutionRange[1]}, "E_rec", "phi_diff");
0122     auto phi_diff_vs_theta_Hist   = momentumDF.Histo2D({"phi_diff_vs_theta", "phi difference vs theta reconstructed; theta reconstructed [rad]; phi difference [deg]", thetaBins, thetaRange[0], thetaRange[1], resolutionBins, phiResolutionRange[0], phiResolutionRange[1]}, "theta_rec", "phi_diff");
0123     auto phi_diff_vs_phi_Hist     = momentumDF.Histo2D({"phi_diff_vs_phi", "phi difference vs phi reconstructed; phi reconstructed [deg]; phi difference [deg]", phiBins, phiRange[0], phiRange[1], resolutionBins, phiResolutionRange[0], phiResolutionRange[1]}, "phi_rec", "phi_diff");
0124 
0125     // Create canvas for momentum component plots
0126     TCanvas *cMomentum = new TCanvas("momentum_canvas", "Momentum Resolution", 3000, 1600);
0127     cMomentum->Divide(3, 2);
0128     cMomentum->cd(1);
0129     px_Hist->Draw("colz");
0130     gPad->SetLogz();
0131     cMomentum->cd(2);
0132     py_Hist->Draw("colz");
0133     gPad->SetLogz();
0134     cMomentum->cd(3);
0135     pz_Hist->Draw("colz");
0136     gPad->SetLogz();
0137     cMomentum->cd(4);
0138     px_diff_Hist->Draw();
0139     cMomentum->cd(5);
0140     py_diff_Hist->Draw();
0141     cMomentum->cd(6);
0142     pz_res_Hist->Draw();
0143     cMomentum->SetGrid();
0144     cMomentum->Update();
0145     // Save the canvas as a PNG file
0146     cMomentum->SaveAs(momentumCanvasName);
0147 
0148     // Create canvas for energy, theta and phi resolution plots
0149     TCanvas *cEnergyThetaPhi = new TCanvas("energy_theta_phi_canvas", "Energy, Theta and Phi Resolution", 3000, 1600);
0150     cEnergyThetaPhi->Divide(3, 2);
0151     cEnergyThetaPhi->cd(1);
0152     E_Hist->Draw("colz");
0153     gPad->SetLogz();
0154     cEnergyThetaPhi->cd(2);
0155     theta_Hist->Draw("colz");
0156     gPad->SetLogz();
0157     cEnergyThetaPhi->cd(3);
0158     phi_Hist->Draw("colz");
0159     gPad->SetLogz();
0160     cEnergyThetaPhi->cd(4);
0161     E_res_Hist->Draw();
0162     cEnergyThetaPhi->cd(5);
0163     theta_diff_Hist->Draw();
0164     cEnergyThetaPhi->cd(6);
0165     phi_diff_Hist->Draw();
0166     cEnergyThetaPhi->SetGrid();
0167     cEnergyThetaPhi->Update();
0168     // Save the canvas as a PNG file
0169     cEnergyThetaPhi->SaveAs(energyThetaPhiCanvasName);
0170     
0171     // Create canvas for resolution vs MC values
0172     TCanvas *cResolutionVsMC = new TCanvas("resolution_vs_mc_canvas", "Resolution vs MC Values", 3000, 1600);
0173     cResolutionVsMC->Divide(3, 3);
0174     cResolutionVsMC->cd(1);
0175     E_res_vs_E_Hist->Draw("colz");
0176     gPad->SetLogz();
0177     cResolutionVsMC->cd(2);
0178     E_res_vs_theta_Hist->Draw("colz");
0179     gPad->SetLogz();
0180     cResolutionVsMC->cd(3);
0181     E_res_vs_phi_Hist->Draw("colz");
0182     gPad->SetLogz();
0183     cResolutionVsMC->cd(4);
0184     theta_diff_vs_E_Hist->Draw("colz");
0185     gPad->SetLogz();
0186     cResolutionVsMC->cd(5);
0187     theta_diff_vs_theta_Hist->Draw("colz");
0188     gPad->SetLogz();
0189     cResolutionVsMC->cd(6);
0190     theta_diff_vs_phi_Hist->Draw("colz");
0191     gPad->SetLogz();
0192     cResolutionVsMC->cd(7);
0193     phi_diff_vs_E_Hist->Draw("colz");
0194     gPad->SetLogz();
0195     cResolutionVsMC->cd(8);
0196     phi_diff_vs_theta_Hist->Draw("colz");
0197     gPad->SetLogz();
0198     cResolutionVsMC->cd(9);
0199     phi_diff_vs_phi_Hist->Draw("colz");
0200     gPad->SetLogz();
0201     cResolutionVsMC->SetGrid();
0202     cResolutionVsMC->Update();
0203     // Save the canvas as a PNG file
0204     cResolutionVsMC->SaveAs(relationCanvasName);
0205 
0206     // Fit Gaussians to the E vs E histogram slices
0207     TObjArray* fitE_vs_E = new TObjArray();
0208     TObjArray* fitTheta_vs_E = new TObjArray();
0209     TObjArray* fitPhi_vs_theta = new TObjArray();
0210     E_res_vs_E_Hist->FitSlicesY(nullptr, 1, -1, 0, "Q", fitE_vs_E);
0211     theta_diff_vs_E_Hist->FitSlicesY(nullptr, 1, -1, 0, "Q", fitTheta_vs_E);
0212     phi_diff_vs_theta_Hist->FitSlicesY(nullptr, 1, -1, 0, "Q", fitPhi_vs_theta);
0213 
0214     // Create graphs containing the fitted means and standard deviations
0215     TH1* hE_vs_E_mean         = (TH1*)fitE_vs_E->At(1);     // mean values (index 1)
0216     TH1* hE_vs_E_stddev       = (TH1*)fitE_vs_E->At(2);     // stddev values (index 2)
0217     TH1* hTheta_vs_E_mean     = (TH1*)fitTheta_vs_E->At(1); // mean values (index 1)
0218     TH1* hTheta_vs_E_stddev   = (TH1*)fitTheta_vs_E->At(2); // stddev values (index 2)
0219     TH1* hPhi_vs_theta_mean   = (TH1*)fitPhi_vs_theta->At(1);   // mean values (index 1)
0220     TH1* hPhi_vs_theta_stddev = (TH1*)fitPhi_vs_theta->At(2);   // stddev values (index 2)
0221 
0222     // Create a canvas for the resolution graphs
0223     TCanvas *cResolutionGraphs = new TCanvas("resolution_graphs_canvas", "Resolution Graphs", 1200, 800);
0224     cResolutionGraphs->Divide(3, 2);
0225     cResolutionGraphs->cd(1);
0226     hE_vs_E_mean->SetTitle("Mean Energy Offset vs E MC; Energy MC [GeV]; Mean Energy Offset [GeV]");
0227     hE_vs_E_mean->SetMarkerStyle(20);
0228     hE_vs_E_mean->SetMarkerColor(kBlue);
0229     hE_vs_E_mean->SetMaximum(0.01); // Adjust maximum for better visibility
0230     hE_vs_E_mean->SetMinimum(-0.01); // Adjust minimum for better visibility
0231     hE_vs_E_mean->Draw();
0232     cResolutionGraphs->cd(4);
0233     hE_vs_E_stddev->SetTitle("Energy Resolution vs E MC; Energy MC [GeV]; Energy Resolution [GeV]");
0234     hE_vs_E_stddev->SetMarkerStyle(20);
0235     hE_vs_E_stddev->SetMarkerColor(kRed);
0236     hE_vs_E_stddev->SetMaximum(0.04); // Adjust maximum for better visibility
0237     hE_vs_E_stddev->SetMinimum(0); // Adjust minimum for better visibility
0238     hE_vs_E_stddev->Draw();
0239     cResolutionGraphs->cd(2);
0240     hTheta_vs_E_mean->SetTitle("Mean Theta Offset vs E MC; Energy MC [GeV]; Mean Theta Offset [rad]");
0241     hTheta_vs_E_mean->SetMarkerStyle(20);
0242     hTheta_vs_E_mean->SetMarkerColor(kBlue);
0243     hTheta_vs_E_mean->SetMaximum(0.0003); // Adjust maximum for better visibility
0244     hTheta_vs_E_mean->SetMinimum(-0.0003); // Adjust minimum for better visibility
0245     hTheta_vs_E_mean->Draw();
0246     cResolutionGraphs->cd(5);
0247     hTheta_vs_E_stddev->SetTitle("Std Dev Theta Offset vs E MC; Energy MC [GeV]; Std Dev Theta Offset [rad]");
0248     hTheta_vs_E_stddev->SetMarkerStyle(20);
0249     hTheta_vs_E_stddev->SetMarkerColor(kRed);
0250     hTheta_vs_E_stddev->SetMaximum(0.0005); // Adjust maximum for better visibility
0251     hTheta_vs_E_stddev->SetMinimum(0); // Adjust minimum for better visibility
0252     hTheta_vs_E_stddev->Draw();
0253     cResolutionGraphs->cd(3);
0254     hPhi_vs_theta_mean->SetTitle("Mean Phi Offset vs theta MC; theta MC [rad]; Mean Phi Offset [rad]");
0255     hPhi_vs_theta_mean->SetMarkerStyle(20);
0256     hPhi_vs_theta_mean->SetMarkerColor(kBlue);
0257     hPhi_vs_theta_mean->SetMaximum(5); // Adjust maximum for better visibility
0258     hPhi_vs_theta_mean->SetMinimum(-5); // Adjust minimum for better visibility
0259     hPhi_vs_theta_mean->Draw();
0260     cResolutionGraphs->cd(6);
0261     hPhi_vs_theta_stddev->SetTitle("Std Dev Phi Offset vs theta MC; theta MC [rad]; Std Dev Phi Offset [rad]");
0262     hPhi_vs_theta_stddev->SetMarkerStyle(20);
0263     hPhi_vs_theta_stddev->SetMarkerColor(kRed);
0264     hPhi_vs_theta_stddev->SetMaximum(60); // Adjust maximum for better visibility
0265     hPhi_vs_theta_stddev->SetMinimum(0); // Adjust minimum for better visibility
0266     hPhi_vs_theta_stddev->Draw();
0267     cResolutionGraphs->SetGrid();
0268     cResolutionGraphs->Update();
0269     // Save the canvas as a PNG file
0270     cResolutionGraphs->SaveAs(resolutionGraphsCanvasName);
0271 
0272     TFile *f = new TFile(outFile,"RECREATE");
0273     cMomentum->Write();
0274     cEnergyThetaPhi->Write();
0275     cResolutionVsMC->Write();
0276     cResolutionGraphs->Write();
0277     px_Hist->Write();
0278     py_Hist->Write();
0279     pz_Hist->Write();
0280     px_diff_Hist->Write();
0281     py_diff_Hist->Write();
0282     pz_res_Hist->Write();
0283     E_Hist->Write();
0284     theta_Hist->Write();
0285     phi_Hist->Write();
0286     E_res_Hist->Write();
0287     theta_diff_Hist->Write();
0288     phi_diff_Hist->Write();
0289     E_res_vs_E_Hist->Write();
0290     E_res_vs_theta_Hist->Write();
0291     E_res_vs_phi_Hist->Write();
0292     theta_diff_vs_E_Hist->Write();
0293     theta_diff_vs_theta_Hist->Write();
0294     theta_diff_vs_phi_Hist->Write();
0295     phi_diff_vs_E_Hist->Write();
0296     phi_diff_vs_theta_Hist->Write();
0297     phi_diff_vs_phi_Hist->Write();
0298 
0299     hE_vs_E_mean->Write();
0300     hE_vs_E_stddev->Write();
0301     hTheta_vs_E_mean->Write();
0302     hTheta_vs_E_stddev->Write();
0303     hPhi_vs_theta_mean->Write();
0304     hPhi_vs_theta_stddev->Write();
0305 
0306     f->Close();
0307 
0308 
0309     
0310     // // Get mean and error on the mean of E, theta and phi resolutions
0311     // double mean_E_res = E_res_Hist->GetMean();
0312     // double mean_theta_res = theta_diff_Hist->GetMean();
0313     // double mean_phi_res = phi_diff_Hist->GetMean();
0314     // double mean_E_res_error = E_res_Hist->GetMeanError();
0315     // double mean_theta_res_error = theta_diff_Hist->GetMeanError();
0316     // double mean_phi_res_error = phi_diff_Hist->GetMeanError();
0317 
0318     // // Get standard deviation of E, theta and phi resolutions
0319     // double stddev_E_res = E_res_Hist->GetStdDev();
0320     // double stddev_theta_res = theta_diff_Hist->GetStdDev();
0321     // double stddev_phi_res = phi_diff_Hist->GetStdDev();
0322     // double stddev_E_res_error = E_res_Hist->GetStdDevError();
0323     // double stddev_theta_res_error = theta_diff_Hist->GetStdDevError();
0324     // double stddev_phi_res_error = phi_diff_Hist->GetStdDevError();
0325 
0326     // // Print the resolutions
0327     // std::cout << "Mean E offset: " << mean_E_res << " +/- " << mean_E_res_error << std::endl;
0328     // std::cout << "Mean theta offset: " << mean_theta_res << " +/- " << mean_theta_res_error << std::endl;
0329     // std::cout << "Mean phi offset: " << mean_phi_res << " +/- " << mean_phi_res_error << std::endl;
0330     // std::cout << "Standard deviation of E resolution: " << stddev_E_res << " +/- " << stddev_E_res_error << std::endl;
0331     // std::cout << "Standard deviation of theta resolution: " << stddev_theta_res << " +/- " << stddev_theta_res_error << std::endl;
0332     // std::cout << "Standard deviation of phi resolution: " << stddev_phi_res << " +/- " << stddev_phi_res_error << std::endl;
0333 
0334     // int pass = 0;
0335 
0336     // // Fail if mean is more than 20% of the standard deviation away from zero
0337     // if(std::abs(mean_E_res) > 0.2 * stddev_E_res) {
0338     //     std::cout << "Mean E offset is more than 20\% (" << 0.2 * stddev_E_res << ") of the standard deviation away from zero!" << std::endl;
0339     //     pass = 1;
0340     // }
0341     // if(std::abs(mean_theta_res) > 0.2 * stddev_theta_res) {
0342     //     std::cout << "Mean theta offset is more than 20\% (" << 0.2 * stddev_theta_res << ") of the standard deviation away from zero!" << std::endl;
0343     //     pass = 1;
0344     // }
0345     // if(std::abs(mean_phi_res) > 0.2 * stddev_phi_res) {
0346     //     std::cout << "Mean phi offset is more than 20\% (" << 0.2 * stddev_phi_res << ") of the standard deviation away from zero!" << std::endl;
0347     //     pass = 1;
0348     // }
0349 
0350     // // Resolution limits
0351     // double E_res_limit = 0.05; // 5% resolution
0352     // double theta_res_limit = 0.0001; // 1 mrad resolution
0353     // double phi_res_limit = 30; // 30 degrees resolution
0354 
0355     // // Fail if standard deviation is more than the limit
0356     // if(std::abs(stddev_E_res) > E_res_limit) {
0357     //     std::cout << "Standard deviation of E resolution is more than the limit of " << E_res_limit << "!" << std::endl;
0358     //     pass = 1;
0359     // }
0360     // if(std::abs(stddev_theta_res) > theta_res_limit) {
0361     //     std::cout << "Standard deviation of theta resolution is more than the limit of " << theta_res_limit << " radians!" << std::endl;
0362     //     pass = 1;
0363     // }
0364     // if(std::abs(stddev_phi_res) > phi_res_limit) {
0365     //     std::cout << "Standard deviation of phi resolution is more than the limit of " << phi_res_limit << " degrees!" << std::endl;
0366     //     pass = 1;
0367     // }
0368 
0369 }
0370
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!