From 7c18a23ee0ac29e16bc17f85cee004bfe2e8753f Mon Sep 17 00:00:00 2001
From: Wouter Deconinck <wdconinc@gmail.com>
Date: Mon, 7 Feb 2022 21:52:02 +0000
Subject: [PATCH] Add Sigma and eSigma inclusive kinematics
---
benchmarks/dis/analysis/dis_electrons.cxx | 76 ++++++++++++++++++++++-
options/reconstruction.py | 14 +++++
2 files changed, 89 insertions(+), 1 deletion(-)
diff --git a/benchmarks/dis/analysis/dis_electrons.cxx b/benchmarks/dis/analysis/dis_electrons.cxx
index d1ddb07d..fed2cafc 100644
--- a/benchmarks/dis/analysis/dis_electrons.cxx
+++ b/benchmarks/dis/analysis/dis_electrons.cxx
@@ -80,16 +80,24 @@ int dis_electrons(const std::string& config_name)
.Define("Q2_el", "InclusiveKinematicsElectron.Q2")
.Define("Q2_jb", "InclusiveKinematicsJB.Q2")
.Define("Q2_da", "InclusiveKinematicsDA.Q2")
+ .Define("Q2_sigma", "InclusiveKinematicsSigma.Q2")
+ .Define("Q2_esigma", "InclusiveKinematicseSigma.Q2")
.Define("Q2_el_res", combinatorial_diff_ratio, {"Q2_sim", "Q2_el"})
.Define("Q2_jb_res", combinatorial_diff_ratio, {"Q2_sim", "Q2_jb"})
.Define("Q2_da_res", combinatorial_diff_ratio, {"Q2_sim", "Q2_da"})
+ .Define("Q2_sigma_res", combinatorial_diff_ratio, {"Q2_sim", "Q2_sigma"})
+ .Define("Q2_esigma_res", combinatorial_diff_ratio, {"Q2_sim", "Q2_esigma"})
.Define("x_sim", "InclusiveKinematicsTruth.x")
.Define("x_el", "InclusiveKinematicsElectron.x")
.Define("x_jb", "InclusiveKinematicsJB.x")
.Define("x_da", "InclusiveKinematicsDA.x")
+ .Define("x_sigma", "InclusiveKinematicsSigma.x")
+ .Define("x_esigma", "InclusiveKinematicseSigma.x")
.Define("x_el_res", combinatorial_diff_ratio, {"x_sim", "x_el"})
.Define("x_jb_res", combinatorial_diff_ratio, {"x_sim", "x_jb"})
.Define("x_da_res", combinatorial_diff_ratio, {"x_sim", "x_da"})
+ .Define("x_sigma_res", combinatorial_diff_ratio, {"x_sim", "x_sigma"})
+ .Define("x_esigma_res", combinatorial_diff_ratio, {"x_sim", "x_esigma"})
;
//Q2
@@ -97,17 +105,25 @@ int dis_electrons(const std::string& config_name)
auto h_Q2_el = d0.Histo1D({"h_Q2_el", "; GeV^2; counts", 100, -5, 25}, "Q2_el");
auto h_Q2_jb = d0.Histo1D({"h_Q2_jb", "; GeV^2; counts", 100, -5, 25}, "Q2_jb");
auto h_Q2_da = d0.Histo1D({"h_Q2_da", "; GeV^2; counts", 100, -5, 25}, "Q2_da");
+ auto h_Q2_sigma = d0.Histo1D({"h_Q2_sigma", "; GeV^2; counts", 100, -5, 25}, "Q2_sigma");
+ auto h_Q2_esigma = d0.Histo1D({"h_Q2_esigma", "; GeV^2; counts", 100, -5, 25}, "Q2_esigma");
auto h_Q2_el_res = d0.Histo1D({"h_Q2_el_res", "; ; counts", 100, -1, 1}, "Q2_el_res");
auto h_Q2_jb_res = d0.Histo1D({"h_Q2_jb_res", "; ; counts", 100, -1, 1}, "Q2_jb_res");
auto h_Q2_da_res = d0.Histo1D({"h_Q2_da_res", "; ; counts", 100, -1, 1}, "Q2_da_res");
+ auto h_Q2_sigma_res = d0.Histo1D({"h_Q2_sigma_res", "; ; counts", 100, -1, 1}, "Q2_sigma_res");
+ auto h_Q2_esigma_res = d0.Histo1D({"h_Q2_esigma_res", "; ; counts", 100, -1, 1}, "Q2_esigma_res");
//x
auto h_x_sim = d0.Histo1D({"h_x_sim", "; ; counts", 100, 0, +1}, "x_sim");
auto h_x_el = d0.Histo1D({"h_x_el", "; ; counts", 100, 0, +1}, "x_el");
auto h_x_jb = d0.Histo1D({"h_x_jb", "; ; counts", 100, 0, +1}, "x_jb");
auto h_x_da = d0.Histo1D({"h_x_da", "; ; counts", 100, 0, +1}, "x_da");
+ auto h_x_sigma = d0.Histo1D({"h_x_sigma", "; ; counts", 100, 0, +1}, "x_sigma");
+ auto h_x_esigma = d0.Histo1D({"h_x_esigma", "; ; counts", 100, 0, +1}, "x_esigma");
auto h_x_el_res = d0.Histo1D({"h_x_el_res", "; ; counts", 100, -1, 1}, "x_el_res");
auto h_x_jb_res = d0.Histo1D({"h_x_jb_res", "; ; counts", 100, -1, 1}, "x_jb_res");
auto h_x_da_res = d0.Histo1D({"h_x_da_res", "; ; counts", 100, -1, 1}, "x_da_res");
+ auto h_x_sigma_res = d0.Histo1D({"h_x_sigma_res", "; ; counts", 100, -1, 1}, "x_sigma_res");
+ auto h_x_esigma_res = d0.Histo1D({"h_x_esigma_res", "; ; counts", 100, -1, 1}, "x_esigma_res");
TFitResultPtr f_Q2_el_res = h_Q2_el_res->Fit("gaus", "S");
if (f_Q2_el_res == 0) f_Q2_el_res->Print("V");
@@ -124,6 +140,16 @@ int dis_electrons(const std::string& config_name)
TFitResultPtr f_x_da_res = h_x_da_res->Fit("gaus", "S");
if (f_x_da_res == 0) f_x_da_res->Print("V");
+ TFitResultPtr f_Q2_sigma_res = h_Q2_sigma_res->Fit("gaus", "S");
+ if (f_Q2_sigma_res == 0) f_Q2_sigma_res->Print("V");
+ TFitResultPtr f_x_sigma_res = h_x_sigma_res->Fit("gaus", "S");
+ if (f_x_sigma_res == 0) f_x_sigma_res->Print("V");
+
+ TFitResultPtr f_Q2_esigma_res = h_Q2_esigma_res->Fit("gaus", "S");
+ if (f_Q2_esigma_res == 0) f_Q2_esigma_res->Print("V");
+ TFitResultPtr f_x_esigma_res = h_x_esigma_res->Fit("gaus", "S");
+ if (f_x_esigma_res == 0) f_x_esigma_res->Print("V");
+
// Plot our histograms.
// TODO: to start I'm explicitly plotting the histograms, but want to
// factorize out the plotting code moving forward.
@@ -138,6 +164,8 @@ int dis_electrons(const std::string& config_name)
auto& h2 = *h_Q2_el;
auto& h3 = *h_Q2_jb;
auto& h4 = *h_Q2_da;
+ auto& h5 = *h_Q2_sigma;
+ auto& h6 = *h_Q2_esigma;
// histogram style
h1.SetLineColor(common_bench::plot::kMpBlue);
h1.SetLineWidth(2);
@@ -147,6 +175,10 @@ int dis_electrons(const std::string& config_name)
h3.SetLineWidth(2);
h4.SetLineColor(common_bench::plot::kMpGreen);
h4.SetLineWidth(2);
+ h5.SetLineColor(common_bench::plot::kMpMoss);
+ h5.SetLineWidth(2);
+ h6.SetLineColor(common_bench::plot::kMpCyan);
+ h6.SetLineWidth(2);
// axes
h1.GetXaxis()->CenterTitle();
h1.GetYaxis()->CenterTitle();
@@ -155,6 +187,8 @@ int dis_electrons(const std::string& config_name)
h2.DrawClone("hist same");
h3.DrawClone("hist same");
h4.DrawClone("hist same");
+ h5.DrawClone("hist same");
+ h6.DrawClone("hist same");
// legend
common_bench::plot::draw_label(ebeam, pbeam, detector);
TText* tptr1;
@@ -164,12 +198,16 @@ int dis_electrons(const std::string& config_name)
t1.SetTextSize(25);
tptr1 = t1.AddText("simulated");
tptr1->SetTextColor(common_bench::plot::kMpBlue);
- tptr1 = t1.AddText("EL method");
+ tptr1 = t1.AddText("e method");
tptr1->SetTextColor(common_bench::plot::kMpOrange);
tptr1 = t1.AddText("JB method");
tptr1->SetTextColor(common_bench::plot::kMpRed);
tptr1 = t1.AddText("DA method");
tptr1->SetTextColor(common_bench::plot::kMpGreen);
+ tptr1 = t1.AddText("#Sigma method");
+ tptr1->SetTextColor(common_bench::plot::kMpMoss);
+ tptr1 = t1.AddText("e#Sigma method");
+ tptr1->SetTextColor(common_bench::plot::kMpCyan);
t1.Draw();
c.Print(fmt::format("{}_Q2.png", output_prefix).c_str());
}
@@ -183,6 +221,8 @@ int dis_electrons(const std::string& config_name)
auto& h1 = *h_Q2_el_res;
auto& h2 = *h_Q2_jb_res;
auto& h3 = *h_Q2_da_res;
+ auto& h4 = *h_Q2_sigma_res;
+ auto& h5 = *h_Q2_esigma_res;
// histogram style
h1.SetLineColor(common_bench::plot::kMpOrange);
h1.SetLineWidth(2);
@@ -190,6 +230,10 @@ int dis_electrons(const std::string& config_name)
h2.SetLineWidth(2);
h3.SetLineColor(common_bench::plot::kMpGreen);
h3.SetLineWidth(2);
+ h4.SetLineColor(common_bench::plot::kMpMoss);
+ h4.SetLineWidth(2);
+ h5.SetLineColor(common_bench::plot::kMpCyan);
+ h5.SetLineWidth(2);
// axes
h1.GetXaxis()->CenterTitle();
h1.GetYaxis()->CenterTitle();
@@ -197,6 +241,8 @@ int dis_electrons(const std::string& config_name)
h1.DrawClone("hist");
h2.DrawClone("hist same");
h3.DrawClone("hist same");
+ h4.DrawClone("hist same");
+ h5.DrawClone("hist same");
// legend
common_bench::plot::draw_label(ebeam, pbeam, detector);
TText* tptr1;
@@ -210,6 +256,10 @@ int dis_electrons(const std::string& config_name)
tptr1->SetTextColor(common_bench::plot::kMpRed);
tptr1 = t1.AddText("DA method");
tptr1->SetTextColor(common_bench::plot::kMpGreen);
+ tptr1 = t1.AddText("#Sigma method");
+ tptr1->SetTextColor(common_bench::plot::kMpMoss);
+ tptr1 = t1.AddText("e#Sigma method");
+ tptr1->SetTextColor(common_bench::plot::kMpCyan);
t1.Draw();
c.Print(fmt::format("{}_Q2_resolution.png", output_prefix).c_str());
}
@@ -224,6 +274,8 @@ int dis_electrons(const std::string& config_name)
auto& h2 = *h_x_el;
auto& h3 = *h_x_jb;
auto& h4 = *h_x_da;
+ auto& h5 = *h_x_sigma;
+ auto& h6 = *h_x_esigma;
// histogram style
h1.SetLineColor(common_bench::plot::kMpBlue);
h1.SetLineWidth(2);
@@ -233,6 +285,10 @@ int dis_electrons(const std::string& config_name)
h3.SetLineWidth(2);
h4.SetLineColor(common_bench::plot::kMpGreen);
h4.SetLineWidth(2);
+ h5.SetLineColor(common_bench::plot::kMpMoss);
+ h5.SetLineWidth(2);
+ h6.SetLineColor(common_bench::plot::kMpCyan);
+ h6.SetLineWidth(2);
// axes
h1.GetXaxis()->CenterTitle();
h1.GetYaxis()->CenterTitle();
@@ -241,6 +297,8 @@ int dis_electrons(const std::string& config_name)
h2.DrawClone("hist same");
h3.DrawClone("hist same");
h4.DrawClone("hist same");
+ h5.DrawClone("hist same");
+ h6.DrawClone("hist same");
// legend
common_bench::plot::draw_label(ebeam, pbeam, detector);
TText* tptr1;
@@ -256,6 +314,10 @@ int dis_electrons(const std::string& config_name)
tptr1->SetTextColor(common_bench::plot::kMpRed);
tptr1 = t1.AddText("DA method");
tptr1->SetTextColor(common_bench::plot::kMpGreen);
+ tptr1 = t1.AddText("#Sigma method");
+ tptr1->SetTextColor(common_bench::plot::kMpMoss);
+ tptr1 = t1.AddText("e#Sigma method");
+ tptr1->SetTextColor(common_bench::plot::kMpCyan);
t1.Draw();
c.Print(fmt::format("{}_x.png", output_prefix).c_str());
}
@@ -269,6 +331,8 @@ int dis_electrons(const std::string& config_name)
auto& h1 = *h_x_el_res;
auto& h2 = *h_x_jb_res;
auto& h3 = *h_x_da_res;
+ auto& h4 = *h_x_sigma_res;
+ auto& h5 = *h_x_esigma_res;
// histogram style
h1.SetLineColor(common_bench::plot::kMpOrange);
h1.SetLineWidth(2);
@@ -276,6 +340,10 @@ int dis_electrons(const std::string& config_name)
h2.SetLineWidth(2);
h3.SetLineColor(common_bench::plot::kMpGreen);
h3.SetLineWidth(2);
+ h4.SetLineColor(common_bench::plot::kMpMoss);
+ h4.SetLineWidth(2);
+ h5.SetLineColor(common_bench::plot::kMpCyan);
+ h5.SetLineWidth(2);
// axes
h1.GetXaxis()->CenterTitle();
h1.GetYaxis()->CenterTitle();
@@ -283,6 +351,8 @@ int dis_electrons(const std::string& config_name)
h1.DrawClone("hist");
h2.DrawClone("hist same");
h3.DrawClone("hist same");
+ h4.DrawClone("hist same");
+ h5.DrawClone("hist same");
// legend
common_bench::plot::draw_label(ebeam, pbeam, detector);
TText* tptr1;
@@ -296,6 +366,10 @@ int dis_electrons(const std::string& config_name)
tptr1->SetTextColor(common_bench::plot::kMpRed);
tptr1 = t1.AddText("DA method");
tptr1->SetTextColor(common_bench::plot::kMpGreen);
+ tptr1 = t1.AddText("#Sigma method");
+ tptr1->SetTextColor(common_bench::plot::kMpMoss);
+ tptr1 = t1.AddText("e#Sigma method");
+ tptr1->SetTextColor(common_bench::plot::kMpCyan);
t1.Draw();
c.Print(fmt::format("{}_x_resolution.png", output_prefix).c_str());
}
diff --git a/options/reconstruction.py b/options/reconstruction.py
index 1b5c6228..178e31f2 100644
--- a/options/reconstruction.py
+++ b/options/reconstruction.py
@@ -123,6 +123,8 @@ from Configurables import Jug__Reco__ParticlesFromTrackFit as ParticlesFromTrack
from Configurables import Jug__Reco__InclusiveKinematicsElectron as InclusiveKinematicsElectron
from Configurables import Jug__Reco__InclusiveKinematicsDA as InclusiveKinematicsDA
from Configurables import Jug__Reco__InclusiveKinematicsJB as InclusiveKinematicsJB
+from Configurables import Jug__Reco__InclusiveKinematicsSigma as InclusiveKinematicsSigma
+from Configurables import Jug__Reco__InclusiveKinematicseSigma as InclusiveKinematicseSigma
from Configurables import Jug__Reco__FarForwardParticles as FFRecoRP
from Configurables import Jug__Reco__FarForwardParticlesOMD as FFRecoOMD
@@ -764,6 +766,18 @@ incl_kin_da = InclusiveKinematicsDA("incl_kin_da",
outputData="InclusiveKinematicsDA"
)
algorithms.append(incl_kin_da)
+incl_kin_sigma = InclusiveKinematicsSigma("incl_kin_sigma",
+ inputMCParticles="mcparticles",
+ inputParticles="ReconstructedParticles",
+ outputData="InclusiveKinematicsSigma"
+)
+algorithms.append(incl_kin_sigma)
+incl_kin_esigma = InclusiveKinematicseSigma("incl_kin_esigma",
+ inputMCParticles="mcparticles",
+ inputParticles="ReconstructedParticles",
+ outputData="InclusiveKinematicseSigma"
+)
+algorithms.append(incl_kin_esigma)
# Output
podout = PodioOutput("out", filename=output_rec)
--
GitLab