Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-01-31 10:30:15 

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