Skip to content
Snippets Groups Projects
Commit 69a7be56 authored by Wouter Deconinck's avatar Wouter Deconinck
Browse files

Inclusive Kinematics JB DA

parent d3b07bcb
No related branches found
No related tags found
1 merge request!113Inclusive Kinematics JB DA
...@@ -73,26 +73,42 @@ int dis_electrons(const std::string& config_name) ...@@ -73,26 +73,42 @@ int dis_electrons(const std::string& config_name)
}; };
auto d0 = d.Define("Q2_sim", "InclusiveKinematicsTruth.Q2") auto d0 = d.Define("Q2_sim", "InclusiveKinematicsTruth.Q2")
.Define("Q2_rec", "InclusiveKinematicsElectron.Q2") .Define("Q2_el", "InclusiveKinematicsElectron.Q2")
.Define("Q2_res", combinatorial_diff_ratio, {"Q2_sim", "Q2_rec"}) .Define("Q2_jb", "InclusiveKinematicsJB.Q2")
.Define("Q2_da", "InclusiveKinematicsDA.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("x_sim", "InclusiveKinematicsTruth.x") .Define("x_sim", "InclusiveKinematicsTruth.x")
.Define("x_rec", "InclusiveKinematicsElectron.x") .Define("x_el", "InclusiveKinematicsElectron.x")
.Define("x_res", combinatorial_diff_ratio, {"x_sim", "x_rec"}) .Define("x_jb", "InclusiveKinematicsJB.x")
.Define("x_da", "InclusiveKinematicsDA.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"})
; ;
//Q2 //Q2
auto h_Q2_sim = d0.Histo1D({"h_Q2_sim", "; GeV^2; counts", 100, -5, 25}, "Q2_sim"); auto h_Q2_sim = d0.Histo1D({"h_Q2_sim", "; GeV^2; counts", 100, -5, 25}, "Q2_sim");
auto h_Q2_rec = d0.Histo1D({"h_Q2_rec", "; GeV^2; counts", 100, -5, 25}, "Q2_rec"); auto h_Q2_el = d0.Histo1D({"h_Q2_el", "; GeV^2; counts", 100, -5, 25}, "Q2_el");
auto h_Q2_res = d0.Histo1D({"h_Q2_res", "; ; counts", 100, -1, 1}, "Q2_res"); 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_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");
//x //x
auto h_x_sim = d0.Histo1D({"h_x_sim", "; ; counts", 100, 0, +1}, "x_sim"); auto h_x_sim = d0.Histo1D({"h_x_sim", "; ; counts", 100, 0, +1}, "x_sim");
auto h_x_rec = d0.Histo1D({"h_x_rec", "; ; counts", 100, 0, +1}, "x_rec"); auto h_x_el = d0.Histo1D({"h_x_el", "; ; counts", 100, 0, +1}, "x_el");
auto h_x_res = d0.Histo1D({"h_x_res", "; ; counts", 100, -1, 1}, "x_res"); 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_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");
TFitResultPtr f_Q2_res = h_Q2_res->Fit("gaus", "S"); TFitResultPtr f_Q2_el_res = h_Q2_el_res->Fit("gaus", "S");
if (f_Q2_res == 0) f_Q2_res->Print("V"); if (f_Q2_el_res == 0) f_Q2_el_res->Print("V");
TFitResultPtr f_x_res = h_x_res->Fit("gaus", "S"); TFitResultPtr f_x_el_res = h_x_el_res->Fit("gaus", "S");
if (f_x_res == 0) f_x_res->Print("V"); if (f_x_el_res == 0) f_x_el_res->Print("V");
// Plot our histograms. // Plot our histograms.
// TODO: to start I'm explicitly plotting the histograms, but want to // TODO: to start I'm explicitly plotting the histograms, but want to
...@@ -105,18 +121,26 @@ int dis_electrons(const std::string& config_name) ...@@ -105,18 +121,26 @@ int dis_electrons(const std::string& config_name)
gPad->SetLogx(false); gPad->SetLogx(false);
gPad->SetLogy(true); gPad->SetLogy(true);
auto& h1 = *h_Q2_sim; auto& h1 = *h_Q2_sim;
auto& h2 = *h_Q2_rec; auto& h2 = *h_Q2_el;
auto& h3 = *h_Q2_jb;
auto& h4 = *h_Q2_da;
// histogram style // histogram style
h1.SetLineColor(common_bench::plot::kMpBlue); h1.SetLineColor(common_bench::plot::kMpBlue);
h1.SetLineWidth(2); h1.SetLineWidth(2);
h2.SetLineColor(common_bench::plot::kMpOrange); h2.SetLineColor(common_bench::plot::kMpOrange);
h2.SetLineWidth(2); h2.SetLineWidth(2);
h3.SetLineColor(common_bench::plot::kMpRed);
h3.SetLineWidth(2);
h4.SetLineColor(common_bench::plot::kMpGreen);
h4.SetLineWidth(2);
// axes // axes
h1.GetXaxis()->CenterTitle(); h1.GetXaxis()->CenterTitle();
h1.GetYaxis()->CenterTitle(); h1.GetYaxis()->CenterTitle();
// draw everything // draw everything
h1.DrawClone("hist"); h1.DrawClone("hist");
h2.DrawClone("hist same"); h2.DrawClone("hist same");
h3.DrawClone("hist same");
h4.DrawClone("hist same");
// FIXME hardcoded beam configuration // FIXME hardcoded beam configuration
common_bench::plot::draw_label(18, 275, detector); common_bench::plot::draw_label(18, 275, detector);
TText* tptr1; TText* tptr1;
...@@ -126,8 +150,12 @@ int dis_electrons(const std::string& config_name) ...@@ -126,8 +150,12 @@ int dis_electrons(const std::string& config_name)
t1.SetTextSize(25); t1.SetTextSize(25);
tptr1 = t1.AddText("simulated"); tptr1 = t1.AddText("simulated");
tptr1->SetTextColor(common_bench::plot::kMpBlue); tptr1->SetTextColor(common_bench::plot::kMpBlue);
tptr1 = t1.AddText("reconstructed"); tptr1 = t1.AddText("electron");
tptr1->SetTextColor(common_bench::plot::kMpOrange); tptr1->SetTextColor(common_bench::plot::kMpOrange);
tptr1 = t1.AddText("JB");
tptr1->SetTextColor(common_bench::plot::kMpRed);
tptr1 = t1.AddText("DA");
tptr1->SetTextColor(common_bench::plot::kMpGreen);
t1.Draw(); t1.Draw();
c.Print(fmt::format("{}Q2.png", output_prefix).c_str()); c.Print(fmt::format("{}Q2.png", output_prefix).c_str());
} }
...@@ -138,15 +166,23 @@ int dis_electrons(const std::string& config_name) ...@@ -138,15 +166,23 @@ int dis_electrons(const std::string& config_name)
c.cd(); c.cd();
gPad->SetLogx(false); gPad->SetLogx(false);
gPad->SetLogy(true); gPad->SetLogy(true);
auto& h1 = *h_Q2_res; auto& h1 = *h_Q2_el_res;
auto& h2 = *h_Q2_jb_res;
auto& h3 = *h_Q2_da_res;
// histogram style // histogram style
h1.SetLineColor(common_bench::plot::kMpBlue); h1.SetLineColor(common_bench::plot::kMpBlue);
h1.SetLineWidth(2); h1.SetLineWidth(2);
h2.SetLineColor(common_bench::plot::kMpOrange);
h2.SetLineWidth(2);
h3.SetLineColor(common_bench::plot::kMpRed);
h3.SetLineWidth(2);
// axes // axes
h1.GetXaxis()->CenterTitle(); h1.GetXaxis()->CenterTitle();
h1.GetYaxis()->CenterTitle(); h1.GetYaxis()->CenterTitle();
// draw everything // draw everything
h1.DrawClone("hist"); h1.DrawClone("hist");
h2.DrawClone("hist same");
h3.DrawClone("hist same");
// FIXME hardcoded beam configuration // FIXME hardcoded beam configuration
common_bench::plot::draw_label(18, 275, detector); common_bench::plot::draw_label(18, 275, detector);
c.Print(fmt::format("{}Q2resolution.png", output_prefix).c_str()); c.Print(fmt::format("{}Q2resolution.png", output_prefix).c_str());
...@@ -159,18 +195,26 @@ int dis_electrons(const std::string& config_name) ...@@ -159,18 +195,26 @@ int dis_electrons(const std::string& config_name)
gPad->SetLogx(true); gPad->SetLogx(true);
gPad->SetLogy(true); gPad->SetLogy(true);
auto& h1 = *h_x_sim; auto& h1 = *h_x_sim;
auto& h2 = *h_x_rec; auto& h2 = *h_x_el;
auto& h3 = *h_x_jb;
auto& h4 = *h_x_da;
// histogram style // histogram style
h1.SetLineColor(common_bench::plot::kMpBlue); h1.SetLineColor(common_bench::plot::kMpBlue);
h1.SetLineWidth(2); h1.SetLineWidth(2);
h2.SetLineColor(common_bench::plot::kMpOrange); h2.SetLineColor(common_bench::plot::kMpOrange);
h2.SetLineWidth(2); h2.SetLineWidth(2);
h3.SetLineColor(common_bench::plot::kMpRed);
h3.SetLineWidth(2);
h4.SetLineColor(common_bench::plot::kMpGreen);
h4.SetLineWidth(2);
// axes // axes
h1.GetXaxis()->CenterTitle(); h1.GetXaxis()->CenterTitle();
h1.GetYaxis()->CenterTitle(); h1.GetYaxis()->CenterTitle();
// draw everything // draw everything
h1.DrawClone("hist"); h1.DrawClone("hist");
h2.DrawClone("hist same"); h2.DrawClone("hist same");
h3.DrawClone("hist same");
h4.DrawClone("hist same");
// FIXME hardcoded beam configuration // FIXME hardcoded beam configuration
common_bench::plot::draw_label(18, 275, detector); common_bench::plot::draw_label(18, 275, detector);
TText* tptr1; TText* tptr1;
...@@ -192,15 +236,23 @@ int dis_electrons(const std::string& config_name) ...@@ -192,15 +236,23 @@ int dis_electrons(const std::string& config_name)
c.cd(); c.cd();
gPad->SetLogx(false); gPad->SetLogx(false);
gPad->SetLogy(true); gPad->SetLogy(true);
auto& h1 = *h_x_res; auto& h1 = *h_x_el_res;
auto& h2 = *h_x_jb_res;
auto& h3 = *h_x_da_res;
// histogram style // histogram style
h1.SetLineColor(common_bench::plot::kMpBlue); h1.SetLineColor(common_bench::plot::kMpBlue);
h1.SetLineWidth(2); h1.SetLineWidth(2);
h2.SetLineColor(common_bench::plot::kMpOrange);
h2.SetLineWidth(2);
h3.SetLineColor(common_bench::plot::kMpRed);
h3.SetLineWidth(2);
// axes // axes
h1.GetXaxis()->CenterTitle(); h1.GetXaxis()->CenterTitle();
h1.GetYaxis()->CenterTitle(); h1.GetYaxis()->CenterTitle();
// draw everything // draw everything
h1.DrawClone("hist"); h1.DrawClone("hist");
h2.DrawClone("hist same");
h3.DrawClone("hist same");
// FIXME hardcoded beam configuration // FIXME hardcoded beam configuration
common_bench::plot::draw_label(18, 275, detector); common_bench::plot::draw_label(18, 275, detector);
c.Print(fmt::format("{}xresolution.png", output_prefix).c_str()); c.Print(fmt::format("{}xresolution.png", output_prefix).c_str());
......
...@@ -93,6 +93,8 @@ from Configurables import Jug__Reco__CKFTracking as CKFTracking ...@@ -93,6 +93,8 @@ from Configurables import Jug__Reco__CKFTracking as CKFTracking
from Configurables import Jug__Reco__ParticlesFromTrackFit as ParticlesFromTrackFit from Configurables import Jug__Reco__ParticlesFromTrackFit as ParticlesFromTrackFit
from Configurables import Jug__Reco__TrajectoryFromTrackFit as TrajectoryFromTrackFit from Configurables import Jug__Reco__TrajectoryFromTrackFit as TrajectoryFromTrackFit
from Configurables import Jug__Reco__InclusiveKinematicsElectron as InclusiveKinematicsElectron 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__FarForwardParticles as FFRecoRP from Configurables import Jug__Reco__FarForwardParticles as FFRecoRP
from Configurables import Jug__Reco__FarForwardParticlesOMD as FFRecoOMD from Configurables import Jug__Reco__FarForwardParticlesOMD as FFRecoOMD
...@@ -679,13 +681,25 @@ if 'acadia' in detector_version: ...@@ -679,13 +681,25 @@ if 'acadia' in detector_version:
outputHitCollection="MRICHRecHits") outputHitCollection="MRICHRecHits")
algorithms.append(mrich_reco) algorithms.append(mrich_reco)
# Electron kinematics # Inclusive kinematics
electron_incl_kin = InclusiveKinematicsElectron("electron_incl_kin", incl_kin_electron = InclusiveKinematicsElectron("incl_kin_electron",
inputMCParticles="mcparticles", inputMCParticles="mcparticles",
inputParticles="ReconstructedParticles", inputParticles="ReconstructedParticles",
outputData="InclusiveKinematicsElectron" outputData="InclusiveKinematicsElectron"
) )
algorithms.append(electron_incl_kin) algorithms.append(incl_kin_electron)
incl_kin_jb = InclusiveKinematicsJB("incl_kin_jb",
inputMCParticles="mcparticles",
inputParticles="ReconstructedParticles",
outputData="InclusiveKinematicsJB"
)
algorithms.append(incl_kin_jb)
incl_kin_da = InclusiveKinematicsDA("incl_kin_da",
inputMCParticles="mcparticles",
inputParticles="ReconstructedParticles",
outputData="InclusiveKinematicsDA"
)
algorithms.append(incl_kin_da)
# Output # Output
podout = PodioOutput("out", filename=output_rec) podout = PodioOutput("out", filename=output_rec)
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment