benchmarks/lowq2_reconstruction/checkResolutions.C

0001 #include "TFile.h"
0002 #include "TH1F.h"
0003 #include "TMath.h"
0004 #include <iostream>
0005 #include <fstream>
0006
0007 int checkResolutions(const TString inputFile="/home/simong/EIC/detector_benchmarks_anl/sim_output/beamline/acceptanceTestcurrent.edm4hep.root", const TString outputFile="test.json") {
0008
0009     int fail = 0;
0010
0011     TFile *file = TFile::Open(inputFile);
0012     if (!file || file->IsZombie()) {
0013         std::cerr << "Error opening file: " << inputFile << std::endl;
0014         return 1; // Return error code
0015     }
0016
0017     TH1F* E_res_Hist = (TH1F*)file->Get("E_res");
0018     TH1F* theta_diff_Hist = (TH1F*)file->Get("theta_diff");
0019     TH1F* phi_diff_Hist = (TH1F*)file->Get("phi_diff");
0020
0021
0022     double mean_E_res = E_res_Hist->GetMean();
0023     double mean_theta_res = theta_diff_Hist->GetMean();
0024     double mean_phi_res = phi_diff_Hist->GetMean();
0025     double mean_E_res_error = E_res_Hist->GetMeanError();
0026     double mean_theta_res_error = theta_diff_Hist->GetMeanError();
0027     double mean_phi_res_error = phi_diff_Hist->GetMeanError();
0028
0029     // Get standard deviation of E, theta and phi resolutions
0030     double stddev_E_res = E_res_Hist->GetStdDev();
0031     double stddev_theta_res = theta_diff_Hist->GetStdDev();
0032     double stddev_phi_res = phi_diff_Hist->GetStdDev();
0033     double stddev_E_res_error = E_res_Hist->GetStdDevError();
0034     double stddev_theta_res_error = theta_diff_Hist->GetStdDevError();
0035     double stddev_phi_res_error = phi_diff_Hist->GetStdDevError();
0036
0037     // Print the resolutions
0038     std::cout << "Mean E offset: " << mean_E_res << " +/- " << mean_E_res_error << std::endl;
0039     std::cout << "Mean theta offset: " << mean_theta_res << " +/- " << mean_theta_res_error << std::endl;
0040     std::cout << "Mean phi offset: " << mean_phi_res << " +/- " << mean_phi_res_error << std::endl;
0041     std::cout << "Standard deviation of E resolution: " << stddev_E_res << " +/- " << stddev_E_res_error << std::endl;
0042     std::cout << "Standard deviation of theta resolution: " << stddev_theta_res << " +/- " << stddev_theta_res_error << std::endl;
0043     std::cout << "Standard deviation of phi resolution: " << stddev_phi_res << " +/- " << stddev_phi_res_error << std::endl;
0044
0045     // Fail if mean is more than 20% of the standard deviation away from zero
0046     if(std::abs(mean_E_res) > 0.2 * stddev_E_res) {
0047         std::cout << "Mean E offset is more than 20\% (" << 0.2 * stddev_E_res << ") of the standard deviation away from zero!" << std::endl;
0048         fail = 1;
0049     }
0050     if(std::abs(mean_theta_res) > 0.2 * stddev_theta_res) {
0051         std::cout << "Mean theta offset is more than 20\% (" << 0.2 * stddev_theta_res << ") of the standard deviation away from zero!" << std::endl;
0052         fail = 1;
0053     }
0054     if(std::abs(mean_phi_res) > 0.2 * stddev_phi_res) {
0055         std::cout << "Mean phi offset is more than 20\% (" << 0.2 * stddev_phi_res << ") of the standard deviation away from zero!" << std::endl;
0056         fail = 1;
0057     }
0058
0059     // Resolution limits
0060     double E_res_limit = 0.05; // 5% resolution
0061     double theta_res_limit = 0.001; // 1 mrad resolution
0062     double phi_res_limit = 30; // 30 degrees resolution
0063
0064     // Fail if standard deviation is more than the limit
0065     if(std::abs(stddev_E_res) > E_res_limit) {
0066         std::cout << "E resolution is more than the limit of " << E_res_limit << "!" << std::endl;
0067         fail = 1;
0068     }
0069     if(std::abs(stddev_theta_res) > theta_res_limit) {
0070         std::cout << "Theta resolution is more than the limit of " << theta_res_limit << " radians!" << std::endl;
0071         fail = 1;
0072     }
0073     if(std::abs(stddev_phi_res) > phi_res_limit) {
0074         std::cout << "Phi resolution is more than the limit of " << phi_res_limit << " degrees!" << std::endl;
0075         fail = 1;
0076     }
0077
0078     // Create json output file containing the resolutions, errors and offsets
0079     std::ofstream jsonFile(outputFile);
0080     if (!jsonFile.is_open()) {
0081         std::cerr << "Error opening output file: " << outputFile << std::endl;
0082         return 1; // Return error code
0083     }
0084     jsonFile << "{\n";
0085     jsonFile << "  \"mean_E_res\": " << mean_E_res << ",\n";
0086     jsonFile << "  \"mean_theta_res\": " << mean_theta_res << ",\n";
0087     jsonFile << "  \"mean_phi_res\": " << mean_phi_res << ",\n";
0088     jsonFile << "  \"mean_E_res_error\": " << mean_E_res_error << ",\n";
0089     jsonFile << "  \"mean_theta_res_error\": " << mean_theta_res_error << ",\n";
0090     jsonFile << "  \"mean_phi_res_error\": " << mean_phi_res_error << ",\n";
0091     jsonFile << "  \"stddev_E_res\": " << stddev_E_res << ",\n";
0092     jsonFile << "  \"stddev_theta_res\": " << stddev_theta_res << ",\n";
0093     jsonFile << "  \"stddev_phi_res\": " << stddev_phi_res << ",\n";
0094     jsonFile << "  \"stddev_E_res_error\": " << stddev_E_res_error << ",\n";
0095     jsonFile << "  \"stddev_theta_res_error\": " << stddev_theta_res_error << ",\n";
0096     jsonFile << "  \"stddev_phi_res_error\": " << stddev_phi_res_error << ",\n";
0097     jsonFile << "  \"fail\": " << fail << "\n";
0098     jsonFile << "}\n";
0099     jsonFile.close();
0100
0101     return 0;
0102
0103 }