Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2026-09-15 09:32:03 

0001 #include "common_bench/benchmark.h"
0002 #include "common_bench/plot.h"
0003 
0004 #include <cmath>
0005 #include <fstream>
0006 #include <iostream>
0007 #include <string>
0008 #include <vector>
0009 #include <algorithm>
0010 #include <utility>
0011 
0012 #include "ROOT/RDataFrame.hxx"
0013 #include <TH1D.h>
0014 #include <TFitResult.h>
0015 #include <TRandom3.h>
0016 #include <TCanvas.h>
0017 
0018 #include "fmt/color.h"
0019 #include "fmt/core.h"
0020 
0021 #include "nlohmann/json.hpp"
0022 #include "edm4eic/InclusiveKinematicsData.h"
0023 #include "edm4eic/ReconstructedParticleData.h"
0024 
0025 int dis_electrons(const std::string& config_name)
0026 {
0027   // read our configuration
0028   std::ifstream  config_file{config_name};
0029   nlohmann::json config;
0030   config_file >> config;
0031 
0032   const std::string rec_file      = config["rec_file"];
0033   const std::string detector      = config["detector"];
0034   const std::string output_prefix = config["output_prefix"];
0035   const std::string test_tag      = config["test_tag"];
0036   const int ebeam                 = config["ebeam"];
0037   const int pbeam                 = config["pbeam"];
0038 
0039   fmt::print(fmt::emphasis::bold | fg(fmt::color::forest_green),
0040              "Running DIS electron analysis...\n");
0041   fmt::print(" - Detector package: {}\n", detector);
0042   fmt::print(" - input file: {}\n", rec_file);
0043   fmt::print(" - output prefix: {}\n", output_prefix);
0044   fmt::print(" - test tag: {}\n", test_tag);
0045   fmt::print(" - ebeam: {}\n", ebeam);
0046   fmt::print(" - pbeam: {}\n", pbeam);
0047 
0048   // create our test definition
0049   // test_tag
0050   common_bench::Test dis_Q2_resolution{
0051       {{"name", fmt::format("{}_Q2_resolution", test_tag)},
0052        {"title", "DIS Q2 resolution"},
0053        {"description",
0054         fmt::format("DIS Q2 resolution with {}, estimated using a Gaussian fit.", detector)},
0055        {"quantity", "resolution (in %)"},
0056        {"target", "0.1"}}};
0057 
0058   // Run this in multi-threaded mode if desired
0059   ROOT::EnableImplicitMT();
0060   // Use ROOT::RDataFrame directly on the events tree/RNTuple
0061   // This works with both TTree (.root) and RNTuple (.rnt.root) formats
0062   auto d = ROOT::RDataFrame("events", rec_file);
0063 
0064   auto combinatorial_diff_ratio = [] (
0065       const ROOT::VecOps::RVec<float>& v1,
0066       const ROOT::VecOps::RVec<float>& v2
0067   ) {
0068     std::vector<float> v;
0069     for (auto& i1: v1) {
0070       for (auto& i2: v2) {
0071         if (i1 != 0) {
0072           v.push_back((i1-i2)/i1);
0073         }
0074       }
0075     }
0076     return v;
0077   };
0078 
0079   auto d0 = d.Define("Q2_sim", "InclusiveKinematicsTruth.Q2")
0080              .Define("Q2_el", "InclusiveKinematicsElectron.Q2")
0081              .Define("Q2_jb", "InclusiveKinematicsJB.Q2")
0082              .Define("Q2_da", "InclusiveKinematicsDA.Q2")
0083              .Define("Q2_sigma", "InclusiveKinematicsSigma.Q2")
0084              .Define("Q2_esigma", "InclusiveKinematicsESigma.Q2")
0085              .Define("logQ2_sim", "log10(Q2_sim)")
0086              .Define("logQ2_el", "log10(Q2_el)")
0087              .Define("logQ2_jb", "log10(Q2_jb)")
0088              .Define("logQ2_da", "log10(Q2_da)")
0089              .Define("logQ2_sigma", "log10(Q2_sigma)")
0090              .Define("logQ2_esigma", "log10(Q2_esigma)")
0091              .Define("Q2_el_res", combinatorial_diff_ratio, {"Q2_sim", "Q2_el"})
0092              .Define("Q2_jb_res", combinatorial_diff_ratio, {"Q2_sim", "Q2_jb"})
0093              .Define("Q2_da_res", combinatorial_diff_ratio, {"Q2_sim", "Q2_da"})
0094              .Define("Q2_sigma_res", combinatorial_diff_ratio, {"Q2_sim", "Q2_sigma"})
0095              .Define("Q2_esigma_res", combinatorial_diff_ratio, {"Q2_sim", "Q2_esigma"})
0096              .Define("x_sim", "InclusiveKinematicsTruth.x")
0097              .Define("x_el", "InclusiveKinematicsElectron.x")
0098              .Define("x_jb", "InclusiveKinematicsJB.x")
0099              .Define("x_da", "InclusiveKinematicsDA.x")
0100              .Define("x_sigma", "InclusiveKinematicsSigma.x")
0101              .Define("x_esigma", "InclusiveKinematicsESigma.x")
0102              .Define("x_el_res", combinatorial_diff_ratio, {"x_sim", "x_el"})
0103              .Define("x_jb_res", combinatorial_diff_ratio, {"x_sim", "x_jb"})
0104              .Define("x_da_res", combinatorial_diff_ratio, {"x_sim", "x_da"})
0105              .Define("x_sigma_res", combinatorial_diff_ratio, {"x_sim", "x_sigma"})
0106              .Define("x_esigma_res", combinatorial_diff_ratio, {"x_sim", "x_esigma"})
0107              ;
0108 
0109   //Q2
0110   auto h_Q2_sim = d0.Histo1D({"h_Q2_sim", "; GeV^2; counts", 100, -5, 25}, "Q2_sim");
0111   auto h_Q2_el = d0.Histo1D({"h_Q2_el", "; GeV^2; counts", 100, -5, 25}, "Q2_el");
0112   auto h_Q2_jb = d0.Histo1D({"h_Q2_jb", "; GeV^2; counts", 100, -5, 25}, "Q2_jb");
0113   auto h_Q2_da = d0.Histo1D({"h_Q2_da", "; GeV^2; counts", 100, -5, 25}, "Q2_da");
0114   auto h_Q2_sigma = d0.Histo1D({"h_Q2_sigma", "; GeV^2; counts", 100, -5, 25}, "Q2_sigma");
0115   auto h_Q2_esigma = d0.Histo1D({"h_Q2_esigma", "; GeV^2; counts", 100, -5, 25}, "Q2_esigma");
0116   auto h_logQ2_sim = d0.Histo1D({"h_logQ2_sim", "; GeV^2; counts", 100, -1, 4}, "logQ2_sim");
0117   auto h_logQ2_el = d0.Histo1D({"h_logQ2_el", "; GeV^2; counts", 100, -1, 4}, "logQ2_el");
0118   auto h_logQ2_jb = d0.Histo1D({"h_logQ2_jb", "; GeV^2; counts", 100, -1, 4}, "logQ2_jb");
0119   auto h_logQ2_da = d0.Histo1D({"h_logQ2_da", "; GeV^2; counts", 100, -1, 4}, "logQ2_da");
0120   auto h_logQ2_sigma = d0.Histo1D({"h_logQ2_sigma", "; GeV^2; counts", 100, -1, 4}, "logQ2_sigma");
0121   auto h_logQ2_esigma = d0.Histo1D({"h_logQ2_esigma", "; GeV^2; counts", 100, -1, 4}, "logQ2_esigma");
0122   auto h_Q2_el_res = d0.Histo1D({"h_Q2_el_res", ";      ; counts", 100, -1,  1}, "Q2_el_res");
0123   auto h_Q2_jb_res = d0.Histo1D({"h_Q2_jb_res", ";      ; counts", 100, -1,  1}, "Q2_jb_res");
0124   auto h_Q2_da_res = d0.Histo1D({"h_Q2_da_res", ";      ; counts", 100, -1,  1}, "Q2_da_res");
0125   auto h_Q2_sigma_res = d0.Histo1D({"h_Q2_sigma_res", ";      ; counts", 100, -1,  1}, "Q2_sigma_res");
0126   auto h_Q2_esigma_res = d0.Histo1D({"h_Q2_esigma_res", ";      ; counts", 100, -1,  1}, "Q2_esigma_res");
0127   //x
0128   auto h_x_sim = d0.Histo1D({"h_x_sim", "; ; counts", 100, 0, +1}, "x_sim");
0129   auto h_x_el = d0.Histo1D({"h_x_el", "; ; counts", 100, 0, +1}, "x_el");
0130   auto h_x_jb = d0.Histo1D({"h_x_jb", "; ; counts", 100, 0, +1}, "x_jb");
0131   auto h_x_da = d0.Histo1D({"h_x_da", "; ; counts", 100, 0, +1}, "x_da");
0132   auto h_x_sigma = d0.Histo1D({"h_x_sigma", "; ; counts", 100, 0, +1}, "x_sigma");
0133   auto h_x_esigma = d0.Histo1D({"h_x_esigma", "; ; counts", 100, 0, +1}, "x_esigma");
0134   auto h_x_el_res = d0.Histo1D({"h_x_el_res", "; ; counts", 100, -1, 1}, "x_el_res");
0135   auto h_x_jb_res = d0.Histo1D({"h_x_jb_res", "; ; counts", 100, -1, 1}, "x_jb_res");
0136   auto h_x_da_res = d0.Histo1D({"h_x_da_res", "; ; counts", 100, -1, 1}, "x_da_res");
0137   auto h_x_sigma_res = d0.Histo1D({"h_x_sigma_res", "; ; counts", 100, -1, 1}, "x_sigma_res");
0138   auto h_x_esigma_res = d0.Histo1D({"h_x_esigma_res", "; ; counts", 100, -1, 1}, "x_esigma_res");
0139 
0140   TFitResultPtr f_Q2_el_res = h_Q2_el_res->Fit("gaus", "S");
0141   if (f_Q2_el_res == 0) f_Q2_el_res->Print("V");
0142   TFitResultPtr f_x_el_res = h_x_el_res->Fit("gaus", "S");
0143   if (f_x_el_res == 0) f_x_el_res->Print("V");
0144 
0145   TFitResultPtr f_Q2_jb_res = h_Q2_jb_res->Fit("gaus", "S");
0146   if (f_Q2_jb_res == 0) f_Q2_jb_res->Print("V");
0147   TFitResultPtr f_x_jb_res = h_x_jb_res->Fit("gaus", "S");
0148   if (f_x_jb_res == 0) f_x_jb_res->Print("V");
0149 
0150   TFitResultPtr f_Q2_da_res = h_Q2_da_res->Fit("gaus", "S");
0151   if (f_Q2_da_res == 0) f_Q2_da_res->Print("V");
0152   TFitResultPtr f_x_da_res = h_x_da_res->Fit("gaus", "S");
0153   if (f_x_da_res == 0) f_x_da_res->Print("V");
0154 
0155   TFitResultPtr f_Q2_sigma_res = h_Q2_sigma_res->Fit("gaus", "S");
0156   if (f_Q2_sigma_res == 0) f_Q2_sigma_res->Print("V");
0157   TFitResultPtr f_x_sigma_res = h_x_sigma_res->Fit("gaus", "S");
0158   if (f_x_sigma_res == 0) f_x_sigma_res->Print("V");
0159 
0160   TFitResultPtr f_Q2_esigma_res = h_Q2_esigma_res->Fit("gaus", "S");
0161   if (f_Q2_esigma_res == 0) f_Q2_esigma_res->Print("V");
0162   TFitResultPtr f_x_esigma_res = h_x_esigma_res->Fit("gaus", "S");
0163   if (f_x_esigma_res == 0) f_x_esigma_res->Print("V");
0164 
0165   // Print summary
0166   fmt::print(fmt::emphasis::bold | fg(fmt::color::forest_green),
0167              "Inclusive kinematics summary:\n");
0168   fmt::print("Q2 resolution:\n");
0169   fmt::print(" - electron mean: {} +/- {}\n",
0170     h_Q2_el_res->GetMean(), h_Q2_el_res->GetMeanError());
0171   fmt::print(" - electron stddev: {} +/- {}\n",
0172     h_Q2_el_res->GetStdDev(), h_Q2_el_res->GetStdDevError());
0173   if (abs(h_Q2_el_res->GetMean()) / h_Q2_el_res->GetMeanError() > 5) {
0174     fmt::print("Q2 electron res mean not 0\n");
0175     return 1;
0176   }
0177   if (f_Q2_el_res == 0) {
0178     fmt::print(" - electron fit: {} +/- {}\n",
0179       f_Q2_el_res->Parameter(1), f_Q2_el_res->Error(1));
0180   } else {
0181     fmt::print("Q2 electron fit failed (FIXME: allowed to fail)\n");
0182     //return 1;
0183   }
0184   fmt::print(" - sigma mean:    {} +/- {}\n",
0185     h_Q2_sigma_res->GetMean(), h_Q2_sigma_res->GetMeanError());
0186   fmt::print(" - sigma stddev:    {} +/- {}\n",
0187     h_Q2_sigma_res->GetStdDev(), h_Q2_sigma_res->GetStdDevError());
0188   if (abs(h_Q2_sigma_res->GetMean()) / h_Q2_sigma_res->GetMeanError() > 5) {
0189     fmt::print("Q2 sigma res mean not 0 (FIXME: allowed to fail)\n");
0190     //return 1;
0191   }
0192   if (f_Q2_sigma_res == 0) {
0193     fmt::print(" - sigma fit:    {} +/- {}\n",
0194       f_Q2_sigma_res->Parameter(1), f_Q2_sigma_res->Error(1));
0195   } else {
0196     fmt::print("Q2 sigma fit failed (FIXME: allowed to fail)\n");
0197     //return 1;
0198   }
0199   fmt::print(" - esigma mean:    {} +/- {}\n",
0200     h_Q2_esigma_res->GetMean(), h_Q2_esigma_res->GetMeanError());
0201   fmt::print(" - esigma stddev:    {} +/- {}\n",
0202     h_Q2_esigma_res->GetStdDev(), h_Q2_esigma_res->GetStdDevError());
0203   if (abs(h_Q2_esigma_res->GetMean()) / h_Q2_esigma_res->GetMeanError() > 5) {
0204     fmt::print("Q2 esigma res mean not 0 (FIXME: allowed to fail)\n");
0205     //return 1;
0206   }
0207   if (f_Q2_esigma_res == 0) {
0208     fmt::print(" - esigma fit:   {} +/- {}\n",
0209       f_Q2_esigma_res->Parameter(1), f_Q2_esigma_res->Error(1));
0210   } else {
0211     fmt::print("Q2 esigma fit failed (FIXME: allowed to fail)\n");
0212     //return 1;
0213   }
0214   fmt::print(" - JB mean:    {} +/- {}\n",
0215     h_Q2_jb_res->GetMean(), h_Q2_jb_res->GetMeanError());
0216   fmt::print(" - JB stddev:    {} +/- {}\n",
0217     h_Q2_jb_res->GetStdDev(), h_Q2_jb_res->GetStdDevError());
0218   if (abs(h_Q2_jb_res->GetMean()) / h_Q2_jb_res->GetMeanError() > 5) {
0219     fmt::print("Q2 JB res mean not 0 (FIXME: allowed to fail)\n");
0220     //return 1;
0221   }
0222   if (f_Q2_jb_res == 0) {
0223     fmt::print(" - JB fit:       {} +/- {}\n",
0224       f_Q2_jb_res->Parameter(1), f_Q2_jb_res->Error(1));
0225   } else {
0226     fmt::print("Q2 JB fit failed (FIXME: allowed to fail)\n");
0227     //return 1;
0228   }
0229   fmt::print(" - DA mean:    {} +/- {}\n",
0230     h_Q2_da_res->GetMean(), h_Q2_da_res->GetMeanError());
0231   fmt::print(" - DA stddev:    {} +/- {}\n",
0232     h_Q2_da_res->GetStdDev(), h_Q2_da_res->GetStdDevError());
0233   if (abs(h_Q2_da_res->GetMean()) / h_Q2_da_res->GetMeanError() > 5) {
0234     fmt::print("Q2 DA res mean not 0 (FIXME: allowed to fail)\n");
0235     //return 1;
0236   }
0237   if (f_Q2_da_res == 0) {
0238     fmt::print(" - DA fit:   (FIXME: allowed to fail)     {} +/- {}\n",
0239       f_Q2_da_res->Parameter(1), f_Q2_da_res->Error(1));
0240   } else {
0241     fmt::print("Q2 DA fit failed (FIXME: allowed to fail)\n");
0242     //return 1;
0243   }
0244   fmt::print("x resolution:\n");
0245   fmt::print(" - electron mean: {} +/- {}\n",
0246     h_x_el_res->GetMean(), h_x_el_res->GetMeanError());
0247   fmt::print(" - electron stddev: {} +/- {}\n",
0248     h_x_el_res->GetStdDev(), h_x_el_res->GetStdDevError());
0249   if (abs(h_x_el_res->GetMean()) / h_x_el_res->GetMeanError() > 5) {
0250     fmt::print("x electron res mean not 0 (FIXME: allowed to fail)\n");
0251     //return 1;
0252   }
0253   if (f_x_el_res == 0) {
0254     fmt::print(" - electron fit: {} +/- {}\n",
0255       f_x_el_res->Parameter(1), f_x_el_res->Error(1));
0256   } else {
0257     fmt::print("x electron fit failed (FIXME: allowed to fail)\n");
0258     //return 1;
0259   }
0260   fmt::print(" - sigma mean: {} +/- {}\n",
0261     h_x_sigma_res->GetMean(), h_x_sigma_res->GetMeanError());
0262   fmt::print(" - sigma stddev: {} +/- {}\n",
0263     h_x_sigma_res->GetStdDev(), h_x_sigma_res->GetStdDevError());
0264   if (abs(h_x_sigma_res->GetMean()) / h_x_sigma_res->GetMeanError() > 5) {
0265     fmt::print("x sigma res mean not 0 (FIXME: allowed to fail)\n");
0266     //return 1;
0267   }
0268   if (f_x_sigma_res == 0) {
0269     fmt::print(" - sigma fit:    {} +/- {}\n",
0270       f_x_sigma_res->Parameter(1), f_x_sigma_res->Error(1));
0271   } else {
0272     fmt::print("x sigma fit failed (FIXME: allowed to fail)\n");
0273     //return 1;
0274   }
0275   fmt::print(" - esigma mean: {} +/- {}\n",
0276     h_x_esigma_res->GetMean(), h_x_esigma_res->GetMeanError());
0277   fmt::print(" - esigma stddev: {} +/- {}\n",
0278     h_x_esigma_res->GetStdDev(), h_x_esigma_res->GetStdDevError());
0279   if (abs(h_x_esigma_res->GetMean()) / h_x_esigma_res->GetMeanError() > 5) {
0280     fmt::print("x esigma res mean not 0 (FIXME: allowed to fail)\n");
0281     //return 1;
0282   }
0283   if (f_x_esigma_res == 0) {
0284     fmt::print(" - esigma fit:   {} +/- {}\n",
0285       f_x_esigma_res->Parameter(1), f_x_esigma_res->Error(1));
0286   } else {
0287     fmt::print("x esigma fit failed (FIXME: allowed to fail)\n");
0288     //return 1;
0289   }
0290   fmt::print(" - JB mean: {} +/- {}\n",
0291     h_x_jb_res->GetMean(), h_x_jb_res->GetMeanError());
0292   fmt::print(" - JB stddev: {} +/- {}\n",
0293     h_x_jb_res->GetStdDev(), h_x_jb_res->GetStdDevError());
0294   if (abs(h_x_jb_res->GetMean()) / h_x_jb_res->GetMeanError() > 5) {
0295     fmt::print("x JB res mean not 0 (FIXME: allowed to fail)\n");
0296     //return 1;
0297   }
0298   if (f_x_jb_res == 0) {
0299     fmt::print(" - JB fit:       {} +/- {}\n",
0300       f_x_jb_res->Parameter(1), f_x_jb_res->Error(1));
0301   } else {
0302     fmt::print("x JB fit failed (FIXME: allowed to fail)\n");
0303     //return 1;
0304   }
0305   fmt::print(" - DA mean: {} +/- {}\n",
0306     h_x_da_res->GetMean(), h_x_da_res->GetMeanError());
0307   fmt::print(" - DA stddev: {} +/- {}\n",
0308     h_x_da_res->GetStdDev(), h_x_da_res->GetStdDevError());
0309   if (abs(h_x_da_res->GetMean()) / h_x_da_res->GetMeanError() > 5) {
0310     fmt::print("x DA res mean not 0 (FIXME: allowed to fail)\n");
0311     //return 1;
0312   }
0313   if (f_x_da_res == 0) {
0314     fmt::print(" - DA fit:       {} +/- {}\n",
0315       f_x_da_res->Parameter(1), f_x_da_res->Error(1));
0316   } else {
0317     fmt::print("x DA fit failed (FIXME: allowed to fail)\n");
0318     //return 1;
0319   }
0320 
0321   // Plot our histograms.
0322   // TODO: to start I'm explicitly plotting the histograms, but want to
0323   // factorize out the plotting code moving forward.
0324 
0325   // Q2 comparison (panels)
0326   {
0327     TCanvas c("c", "c", 1800, 1200);
0328     c.Divide(3,2);
0329     c.cd();
0330     gPad->SetLogx(false);
0331     gPad->SetLogy(true);
0332     auto& h1 = *h_logQ2_sim;
0333     auto& h2 = *h_logQ2_el;
0334     auto& h3 = *h_logQ2_jb;
0335     auto& h4 = *h_logQ2_da;
0336     auto& h5 = *h_logQ2_sigma;
0337     auto& h6 = *h_logQ2_esigma;
0338     // histogram style
0339     h1.SetLineColor(common_bench::plot::kMpBlue);
0340     h1.SetLineWidth(2);
0341     h2.SetLineColor(common_bench::plot::kMpOrange);
0342     h2.SetLineWidth(2);
0343     h3.SetLineColor(common_bench::plot::kMpRed);
0344     h3.SetLineWidth(2);
0345     h4.SetLineColor(common_bench::plot::kMpGreen);
0346     h4.SetLineWidth(2);
0347     h5.SetLineColor(common_bench::plot::kMpMoss);
0348     h5.SetLineWidth(2);
0349     h6.SetLineColor(common_bench::plot::kMpCyan);
0350     h6.SetLineWidth(2);
0351     // axes
0352     h1.GetXaxis()->CenterTitle();
0353     h1.GetYaxis()->CenterTitle();
0354     // draw everything
0355     c.cd(1);
0356     h1.DrawClone("hist");
0357     c.cd(2);
0358     h2.DrawClone("hist");
0359     c.cd(3);
0360     h3.DrawClone("hist");
0361     c.cd(4);
0362     h4.DrawClone("hist");
0363     c.cd(5);
0364     h5.DrawClone("hist");
0365     c.cd(6);
0366     h6.DrawClone("hist");
0367     // legend
0368     common_bench::plot::draw_label(ebeam, pbeam, detector);
0369     TText* tptr1;
0370     TPaveText t1(.6, .8417, .9, .925, "NB NDC");
0371     t1.SetFillColorAlpha(kWhite, 0);
0372     t1.SetTextFont(43);
0373     t1.SetTextSize(25);
0374     tptr1 = t1.AddText("simulated");
0375     tptr1->SetTextColor(common_bench::plot::kMpBlue);
0376     tptr1 = t1.AddText("e method");
0377     tptr1->SetTextColor(common_bench::plot::kMpOrange);
0378     tptr1 = t1.AddText("JB method");
0379     tptr1->SetTextColor(common_bench::plot::kMpRed);
0380     tptr1 = t1.AddText("DA method");
0381     tptr1->SetTextColor(common_bench::plot::kMpGreen);
0382     tptr1 = t1.AddText("#Sigma method");
0383     tptr1->SetTextColor(common_bench::plot::kMpMoss);
0384     tptr1 = t1.AddText("e#Sigma method");
0385     tptr1->SetTextColor(common_bench::plot::kMpCyan);
0386     t1.Draw();
0387     c.Print(fmt::format("{}_logQ2_panels.png", output_prefix).c_str());
0388   }
0389 
0390   // Q2 comparison (overlays)
0391   {
0392     TCanvas c("c", "c", 1200, 1200);
0393     c.cd();
0394     gPad->SetLogx(false);
0395     gPad->SetLogy(true);
0396     auto& h1 = *h_logQ2_sim;
0397     auto& h2 = *h_logQ2_el;
0398     auto& h3 = *h_logQ2_jb;
0399     auto& h4 = *h_logQ2_da;
0400     auto& h5 = *h_logQ2_sigma;
0401     auto& h6 = *h_logQ2_esigma;
0402     // histogram style
0403     h1.SetLineColor(common_bench::plot::kMpBlue);
0404     h1.SetLineWidth(2);
0405     h2.SetLineColor(common_bench::plot::kMpOrange);
0406     h2.SetLineWidth(2);
0407     h3.SetLineColor(common_bench::plot::kMpRed);
0408     h3.SetLineWidth(2);
0409     h4.SetLineColor(common_bench::plot::kMpGreen);
0410     h4.SetLineWidth(2);
0411     h5.SetLineColor(common_bench::plot::kMpMoss);
0412     h5.SetLineWidth(2);
0413     h6.SetLineColor(common_bench::plot::kMpCyan);
0414     h6.SetLineWidth(2);
0415     // axes
0416     h1.GetXaxis()->CenterTitle();
0417     h1.GetYaxis()->CenterTitle();
0418     // draw everything
0419     h1.DrawClone("hist");
0420     h2.DrawClone("hist same");
0421     h3.DrawClone("hist same");
0422     h4.DrawClone("hist same");
0423     h5.DrawClone("hist same");
0424     h6.DrawClone("hist same");
0425     // legend
0426     common_bench::plot::draw_label(ebeam, pbeam, detector);
0427     TText* tptr1;
0428     TPaveText t1(.6, .8417, .9, .925, "NB NDC");
0429     t1.SetFillColorAlpha(kWhite, 0);
0430     t1.SetTextFont(43);
0431     t1.SetTextSize(25);
0432     tptr1 = t1.AddText("simulated");
0433     tptr1->SetTextColor(common_bench::plot::kMpBlue);
0434     tptr1 = t1.AddText("e method");
0435     tptr1->SetTextColor(common_bench::plot::kMpOrange);
0436     tptr1 = t1.AddText("JB method");
0437     tptr1->SetTextColor(common_bench::plot::kMpRed);
0438     tptr1 = t1.AddText("DA method");
0439     tptr1->SetTextColor(common_bench::plot::kMpGreen);
0440     tptr1 = t1.AddText("#Sigma method");
0441     tptr1->SetTextColor(common_bench::plot::kMpMoss);
0442     tptr1 = t1.AddText("e#Sigma method");
0443     tptr1->SetTextColor(common_bench::plot::kMpCyan);
0444     t1.Draw();
0445     c.Print(fmt::format("{}_logQ2_overlays.png", output_prefix).c_str());
0446   }
0447 
0448   // Q2 resolution (panels)
0449   {
0450     TCanvas c("c", "c", 1800, 1200);
0451     c.Divide(3,2);
0452     c.cd();
0453     gPad->SetLogx(false);
0454     gPad->SetLogy(true);
0455     auto& h1 = *h_Q2_el_res;
0456     auto& h2 = *h_Q2_jb_res;
0457     auto& h3 = *h_Q2_da_res;
0458     auto& h4 = *h_Q2_sigma_res;
0459     auto& h5 = *h_Q2_esigma_res;
0460     // histogram style
0461     h1.SetLineColor(common_bench::plot::kMpOrange);
0462     h1.SetLineWidth(2);
0463     h2.SetLineColor(common_bench::plot::kMpRed);
0464     h2.SetLineWidth(2);
0465     h3.SetLineColor(common_bench::plot::kMpGreen);
0466     h3.SetLineWidth(2);
0467     h4.SetLineColor(common_bench::plot::kMpMoss);
0468     h4.SetLineWidth(2);
0469     h5.SetLineColor(common_bench::plot::kMpCyan);
0470     h5.SetLineWidth(2);
0471     // axes
0472     h1.GetXaxis()->CenterTitle();
0473     h1.GetYaxis()->CenterTitle();
0474     // draw everything
0475     c.cd(1);
0476     h1.DrawClone("hist");
0477     c.cd(2);
0478     h2.DrawClone("hist");
0479     c.cd(3);
0480     h3.DrawClone("hist");
0481     c.cd(4);
0482     h4.DrawClone("hist");
0483     c.cd(5);
0484     h5.DrawClone("hist");
0485     // legend
0486     common_bench::plot::draw_label(ebeam, pbeam, detector);
0487     TText* tptr1;
0488     TPaveText t1(.6, .8417, .9, .925, "NB NDC");
0489     t1.SetFillColorAlpha(kWhite, 0);
0490     t1.SetTextFont(43);
0491     t1.SetTextSize(25);
0492     tptr1 = t1.AddText("EL method");
0493     tptr1->SetTextColor(common_bench::plot::kMpOrange);
0494     tptr1 = t1.AddText("JB method");
0495     tptr1->SetTextColor(common_bench::plot::kMpRed);
0496     tptr1 = t1.AddText("DA method");
0497     tptr1->SetTextColor(common_bench::plot::kMpGreen);
0498     tptr1 = t1.AddText("#Sigma method");
0499     tptr1->SetTextColor(common_bench::plot::kMpMoss);
0500     tptr1 = t1.AddText("e#Sigma method");
0501     tptr1->SetTextColor(common_bench::plot::kMpCyan);
0502     t1.Draw();
0503     c.Print(fmt::format("{}_Q2_res_panels.png", output_prefix).c_str());
0504   }
0505 
0506 
0507   // Q2 resolution (overlays)
0508   {
0509     TCanvas c("c", "c", 1200, 1200);
0510     c.cd();
0511     gPad->SetLogx(false);
0512     gPad->SetLogy(true);
0513     auto& h1 = *h_Q2_el_res;
0514     auto& h2 = *h_Q2_jb_res;
0515     auto& h3 = *h_Q2_da_res;
0516     auto& h4 = *h_Q2_sigma_res;
0517     auto& h5 = *h_Q2_esigma_res;
0518     // histogram style
0519     h1.SetLineColor(common_bench::plot::kMpOrange);
0520     h1.SetLineWidth(2);
0521     h2.SetLineColor(common_bench::plot::kMpRed);
0522     h2.SetLineWidth(2);
0523     h3.SetLineColor(common_bench::plot::kMpGreen);
0524     h3.SetLineWidth(2);
0525     h4.SetLineColor(common_bench::plot::kMpMoss);
0526     h4.SetLineWidth(2);
0527     h5.SetLineColor(common_bench::plot::kMpCyan);
0528     h5.SetLineWidth(2);
0529     // axes
0530     h1.GetXaxis()->CenterTitle();
0531     h1.GetYaxis()->CenterTitle();
0532     // draw everything
0533     h1.DrawClone("hist");
0534     h2.DrawClone("hist same");
0535     h3.DrawClone("hist same");
0536     h4.DrawClone("hist same");
0537     h5.DrawClone("hist same");
0538     // legend
0539     common_bench::plot::draw_label(ebeam, pbeam, detector);
0540     TText* tptr1;
0541     TPaveText t1(.6, .8417, .9, .925, "NB NDC");
0542     t1.SetFillColorAlpha(kWhite, 0);
0543     t1.SetTextFont(43);
0544     t1.SetTextSize(25);
0545     tptr1 = t1.AddText("EL method");
0546     tptr1->SetTextColor(common_bench::plot::kMpOrange);
0547     tptr1 = t1.AddText("JB method");
0548     tptr1->SetTextColor(common_bench::plot::kMpRed);
0549     tptr1 = t1.AddText("DA method");
0550     tptr1->SetTextColor(common_bench::plot::kMpGreen);
0551     tptr1 = t1.AddText("#Sigma method");
0552     tptr1->SetTextColor(common_bench::plot::kMpMoss);
0553     tptr1 = t1.AddText("e#Sigma method");
0554     tptr1->SetTextColor(common_bench::plot::kMpCyan);
0555     t1.Draw();
0556     c.Print(fmt::format("{}_Q2_res_overlays.png", output_prefix).c_str());
0557   }
0558 
0559   // x comparison (panels)
0560   {
0561     TCanvas c("c", "c", 1800, 1200);
0562     c.Divide(3,2);
0563     c.cd();
0564     gPad->SetLogx(true);
0565     gPad->SetLogy(true);
0566     auto& h1 = *h_x_sim;
0567     auto& h2 = *h_x_el;
0568     auto& h3 = *h_x_jb;
0569     auto& h4 = *h_x_da;
0570     auto& h5 = *h_x_sigma;
0571     auto& h6 = *h_x_esigma;
0572     // histogram style
0573     h1.SetLineColor(common_bench::plot::kMpBlue);
0574     h1.SetLineWidth(2);
0575     h2.SetLineColor(common_bench::plot::kMpOrange);
0576     h2.SetLineWidth(2);
0577     h3.SetLineColor(common_bench::plot::kMpRed);
0578     h3.SetLineWidth(2);
0579     h4.SetLineColor(common_bench::plot::kMpGreen);
0580     h4.SetLineWidth(2);
0581     h5.SetLineColor(common_bench::plot::kMpMoss);
0582     h5.SetLineWidth(2);
0583     h6.SetLineColor(common_bench::plot::kMpCyan);
0584     h6.SetLineWidth(2);
0585     // axes
0586     h1.GetXaxis()->CenterTitle();
0587     h1.GetYaxis()->CenterTitle();
0588     // draw everything
0589     c.cd(1);
0590     h1.DrawClone("hist");
0591     c.cd(2);
0592     h2.DrawClone("hist");
0593     c.cd(3);
0594     h3.DrawClone("hist");
0595     c.cd(4);
0596     h4.DrawClone("hist");
0597     c.cd(5);
0598     h5.DrawClone("hist");
0599     c.cd(6);
0600     h6.DrawClone("hist");
0601     // legend
0602     common_bench::plot::draw_label(ebeam, pbeam, detector);
0603     TText* tptr1;
0604     TPaveText t1(.6, .8417, .9, .925, "NB NDC");
0605     t1.SetFillColorAlpha(kWhite, 0);
0606     t1.SetTextFont(43);
0607     t1.SetTextSize(25);
0608     tptr1 = t1.AddText("simulated");
0609     tptr1->SetTextColor(common_bench::plot::kMpBlue);
0610     tptr1 = t1.AddText("EL method");
0611     tptr1->SetTextColor(common_bench::plot::kMpOrange);
0612     tptr1 = t1.AddText("JB method");
0613     tptr1->SetTextColor(common_bench::plot::kMpRed);
0614     tptr1 = t1.AddText("DA method");
0615     tptr1->SetTextColor(common_bench::plot::kMpGreen);
0616     tptr1 = t1.AddText("#Sigma method");
0617     tptr1->SetTextColor(common_bench::plot::kMpMoss);
0618     tptr1 = t1.AddText("e#Sigma method");
0619     tptr1->SetTextColor(common_bench::plot::kMpCyan);
0620     t1.Draw();
0621     c.Print(fmt::format("{}_x_panels.png", output_prefix).c_str());
0622   }
0623 
0624   // x comparison (overlays)
0625   {
0626     TCanvas c("c", "c", 1200, 1200);
0627     c.cd();
0628     gPad->SetLogx(true);
0629     gPad->SetLogy(true);
0630     auto& h1 = *h_x_sim;
0631     auto& h2 = *h_x_el;
0632     auto& h3 = *h_x_jb;
0633     auto& h4 = *h_x_da;
0634     auto& h5 = *h_x_sigma;
0635     auto& h6 = *h_x_esigma;
0636     // histogram style
0637     h1.SetLineColor(common_bench::plot::kMpBlue);
0638     h1.SetLineWidth(2);
0639     h2.SetLineColor(common_bench::plot::kMpOrange);
0640     h2.SetLineWidth(2);
0641     h3.SetLineColor(common_bench::plot::kMpRed);
0642     h3.SetLineWidth(2);
0643     h4.SetLineColor(common_bench::plot::kMpGreen);
0644     h4.SetLineWidth(2);
0645     h5.SetLineColor(common_bench::plot::kMpMoss);
0646     h5.SetLineWidth(2);
0647     h6.SetLineColor(common_bench::plot::kMpCyan);
0648     h6.SetLineWidth(2);
0649     // axes
0650     h1.GetXaxis()->CenterTitle();
0651     h1.GetYaxis()->CenterTitle();
0652     // draw everything
0653     h1.DrawClone("hist");
0654     h2.DrawClone("hist same");
0655     h3.DrawClone("hist same");
0656     h4.DrawClone("hist same");
0657     h5.DrawClone("hist same");
0658     h6.DrawClone("hist same");
0659     // legend
0660     common_bench::plot::draw_label(ebeam, pbeam, detector);
0661     TText* tptr1;
0662     TPaveText t1(.6, .8417, .9, .925, "NB NDC");
0663     t1.SetFillColorAlpha(kWhite, 0);
0664     t1.SetTextFont(43);
0665     t1.SetTextSize(25);
0666     tptr1 = t1.AddText("simulated");
0667     tptr1->SetTextColor(common_bench::plot::kMpBlue);
0668     tptr1 = t1.AddText("EL method");
0669     tptr1->SetTextColor(common_bench::plot::kMpOrange);
0670     tptr1 = t1.AddText("JB method");
0671     tptr1->SetTextColor(common_bench::plot::kMpRed);
0672     tptr1 = t1.AddText("DA method");
0673     tptr1->SetTextColor(common_bench::plot::kMpGreen);
0674     tptr1 = t1.AddText("#Sigma method");
0675     tptr1->SetTextColor(common_bench::plot::kMpMoss);
0676     tptr1 = t1.AddText("e#Sigma method");
0677     tptr1->SetTextColor(common_bench::plot::kMpCyan);
0678     t1.Draw();
0679     c.Print(fmt::format("{}_x_overlays.png", output_prefix).c_str());
0680   }
0681 
0682   // x resolution (panels)
0683   {
0684     TCanvas c("c", "c", 1800, 1200);
0685     c.Divide(3,2);
0686     c.cd();
0687     gPad->SetLogx(false);
0688     gPad->SetLogy(true);
0689     auto& h1 = *h_x_el_res;
0690     auto& h2 = *h_x_jb_res;
0691     auto& h3 = *h_x_da_res;
0692     auto& h4 = *h_x_sigma_res;
0693     auto& h5 = *h_x_esigma_res;
0694     // histogram style
0695     h1.SetLineColor(common_bench::plot::kMpOrange);
0696     h1.SetLineWidth(2);
0697     h2.SetLineColor(common_bench::plot::kMpRed);
0698     h2.SetLineWidth(2);
0699     h3.SetLineColor(common_bench::plot::kMpGreen);
0700     h3.SetLineWidth(2);
0701     h4.SetLineColor(common_bench::plot::kMpMoss);
0702     h4.SetLineWidth(2);
0703     h5.SetLineColor(common_bench::plot::kMpCyan);
0704     h5.SetLineWidth(2);
0705     // axes
0706     h1.GetXaxis()->CenterTitle();
0707     h1.GetYaxis()->CenterTitle();
0708     // draw everything
0709     c.cd(1);
0710     h1.DrawClone("hist");
0711     c.cd(2);
0712     h2.DrawClone("hist");
0713     c.cd(3);
0714     h3.DrawClone("hist");
0715     c.cd(4);
0716     h4.DrawClone("hist");
0717     c.cd(5);
0718     h5.DrawClone("hist");
0719     // legend
0720     common_bench::plot::draw_label(ebeam, pbeam, detector);
0721     TText* tptr1;
0722     TPaveText t1(.6, .8417, .9, .925, "NB NDC");
0723     t1.SetFillColorAlpha(kWhite, 0);
0724     t1.SetTextFont(43);
0725     t1.SetTextSize(25);
0726     tptr1 = t1.AddText("EL method");
0727     tptr1->SetTextColor(common_bench::plot::kMpOrange);
0728     tptr1 = t1.AddText("JB method");
0729     tptr1->SetTextColor(common_bench::plot::kMpRed);
0730     tptr1 = t1.AddText("DA method");
0731     tptr1->SetTextColor(common_bench::plot::kMpGreen);
0732     tptr1 = t1.AddText("#Sigma method");
0733     tptr1->SetTextColor(common_bench::plot::kMpMoss);
0734     tptr1 = t1.AddText("e#Sigma method");
0735     tptr1->SetTextColor(common_bench::plot::kMpCyan);
0736     t1.Draw();
0737     c.Print(fmt::format("{}_x_res_panels.png", output_prefix).c_str());
0738   }
0739 
0740   // x resolution (overlays)
0741   {
0742     TCanvas c("c", "c", 1200, 1200);
0743     c.cd();
0744     gPad->SetLogx(false);
0745     gPad->SetLogy(true);
0746     auto& h1 = *h_x_el_res;
0747     auto& h2 = *h_x_jb_res;
0748     auto& h3 = *h_x_da_res;
0749     auto& h4 = *h_x_sigma_res;
0750     auto& h5 = *h_x_esigma_res;
0751     // histogram style
0752     h1.SetLineColor(common_bench::plot::kMpOrange);
0753     h1.SetLineWidth(2);
0754     h2.SetLineColor(common_bench::plot::kMpRed);
0755     h2.SetLineWidth(2);
0756     h3.SetLineColor(common_bench::plot::kMpGreen);
0757     h3.SetLineWidth(2);
0758     h4.SetLineColor(common_bench::plot::kMpMoss);
0759     h4.SetLineWidth(2);
0760     h5.SetLineColor(common_bench::plot::kMpCyan);
0761     h5.SetLineWidth(2);
0762     // axes
0763     h1.GetXaxis()->CenterTitle();
0764     h1.GetYaxis()->CenterTitle();
0765     // draw everything
0766     h1.DrawClone("hist");
0767     h2.DrawClone("hist same");
0768     h3.DrawClone("hist same");
0769     h4.DrawClone("hist same");
0770     h5.DrawClone("hist same");
0771     // legend
0772     common_bench::plot::draw_label(ebeam, pbeam, detector);
0773     TText* tptr1;
0774     TPaveText t1(.6, .8417, .9, .925, "NB NDC");
0775     t1.SetFillColorAlpha(kWhite, 0);
0776     t1.SetTextFont(43);
0777     t1.SetTextSize(25);
0778     tptr1 = t1.AddText("EL method");
0779     tptr1->SetTextColor(common_bench::plot::kMpOrange);
0780     tptr1 = t1.AddText("JB method");
0781     tptr1->SetTextColor(common_bench::plot::kMpRed);
0782     tptr1 = t1.AddText("DA method");
0783     tptr1->SetTextColor(common_bench::plot::kMpGreen);
0784     tptr1 = t1.AddText("#Sigma method");
0785     tptr1->SetTextColor(common_bench::plot::kMpMoss);
0786     tptr1 = t1.AddText("e#Sigma method");
0787     tptr1->SetTextColor(common_bench::plot::kMpCyan);
0788     t1.Draw();
0789     c.Print(fmt::format("{}_x_res_overlays.png", output_prefix).c_str());
0790   }
0791 
0792   common_bench::write_test({dis_Q2_resolution}, fmt::format("{}dis_electrons.json", output_prefix));
0793 
0794   return 0;
0795 }
[ 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!