From 3a0ef3579dd84b33e3cb4bc3a3fc27726c39ad6e Mon Sep 17 00:00:00 2001
From: Ziyue Zhang <zzhan70@uic.edu>
Date: Fri, 22 Jan 2021 16:09:54 -0600
Subject: [PATCH] Add nu Q2 x plots for simulation
---
dvmp/analysis/util.h | 37 +++++++
dvmp/analysis/vm_invar.cxx | 200 +++++++++++++++++++++++++++++++++++++
dvmp/analysis/vm_mass.cxx | 11 +-
dvmp/dvmp.sh | 6 +-
4 files changed, 249 insertions(+), 5 deletions(-)
create mode 100644 dvmp/analysis/vm_invar.cxx
diff --git a/dvmp/analysis/util.h b/dvmp/analysis/util.h
index 0271b787..dad7048e 100644
--- a/dvmp/analysis/util.h
+++ b/dvmp/analysis/util.h
@@ -156,6 +156,43 @@ get_y(const std::pair<ROOT::Math::PxPyPzMVector, ROOT::Math::PxPyPzMVector>&
return 0.5*log((energy_pair + pz_pair)/(energy_pair - pz_pair));
}
+//========================================================================================================
+//for structure functions
+
+struct inv_quant{ //add more when needed
+ double nu, Q2, x;
+};
+
+//for simu
+inline inv_quant calc_inv_quant_simu(const std::vector<ROOT::Math::PxPyPzMVector>& parts){
+ ROOT::Math::PxPyPzMVector q(parts[0] - parts[2]);
+ ROOT::Math::PxPyPzMVector P(parts[3]);
+
+ double nu = q.Dot(P) / P.mass();
+ double Q2 = - q.Dot(q);
+ inv_quant quantities = {nu, Q2, Q2/2./P.mass()/nu};
+ return quantities;
+}
+
+inline double get_nu_simu(inv_quant quantities) {
+ return quantities.nu/1000.;
+}
+inline double get_Q2_simu(inv_quant quantities) {
+ return quantities.Q2;
+}
+inline double get_x_simu(inv_quant quantities) {
+ return quantities.x;
+}
+
+//for tracking, add later
+
+//=========================================================================================================
+
+
+
+
+
+
} // namespace util
#endif
diff --git a/dvmp/analysis/vm_invar.cxx b/dvmp/analysis/vm_invar.cxx
new file mode 100644
index 00000000..f8723733
--- /dev/null
+++ b/dvmp/analysis/vm_invar.cxx
@@ -0,0 +1,200 @@
+#include "benchmark.hh"
+#include "mt.h"
+#include "plot.h"
+#include "util.h"
+
+#include <ROOT/RDataFrame.hxx>
+#include <cmath>
+#include <fmt/color.h>
+#include <fmt/core.h>
+#include <fstream>
+#include <iostream>
+#include <nlohmann/json.hpp>
+#include <string>
+#include <vector>
+
+// Run VM invariant-mass-based benchmarks on an input reconstruction file for
+// a desired vector meson (e.g. jpsi) and a desired decay particle (e.g. muon)
+// Output figures are written to our output prefix (which includes the output
+// file prefix), and labeled with our detector name.
+// TODO: I think it would be better to pass small json configuration file to
+// the test, instead of this ever-expanding list of function arguments.
+// FIXME: MC does not trace back into particle history. Need to fix that
+int vm_invar(const std::string& config_name) {
+ // read our configuration
+ std::ifstream config_file{config_name};
+ nlohmann::json config;
+ config_file >> config;
+
+ const std::string rec_file = config["rec_file"];
+ const std::string vm_name = config["vm_name"];
+ const std::string decay_name = config["decay"];
+ const std::string detector = config["detector"];
+ std::string output_prefix = config["output_prefix"];
+ const std::string test_tag = config["test_tag"];
+
+ fmt::print(fmt::emphasis::bold | fg(fmt::color::forest_green),
+ "Running VM invariant mass analysis...\n");
+ fmt::print(" - Vector meson: {}\n", vm_name);
+ fmt::print(" - Decay particle: {}\n", decay_name);
+ fmt::print(" - Detector package: {}\n", detector);
+ fmt::print(" - output prefix: {}\n", output_prefix);
+
+ // create our test definition
+ // test_tag
+ eic::util::Test vm_mass_resolution_test{
+ {{"name",
+ fmt::format("{}_{}_{}_mass_resolution", test_tag, vm_name, decay_name)},
+ {"title",
+ fmt::format("{} -> {} Invariant Mass Resolution", vm_name, decay_name)},
+ {"description", "Invariant Mass Resolution calculated from raw "
+ "tracking data using a Gaussian fit."},
+ {"quantity", "resolution"},
+ {"target", ".1"}}};
+
+ // Run this in multi-threaded mode if desired
+ ROOT::EnableImplicitMT(kNumThreads);
+
+ // The particles we are looking for. E.g. J/psi decaying into e+e-
+ const double vm_mass = util::get_pdg_mass(vm_name);
+ const double decay_mass = util::get_pdg_mass(decay_name);
+
+ // Ensure our output prefix always ends on a dot, a slash or a dash
+ if (output_prefix.back() != '.' && output_prefix.back() != '/' &&
+ output_prefix.back() != '-') {
+ output_prefix += "-";
+ }
+
+ // Open our input file file as a dataframe
+ ROOT::RDataFrame d{"events", rec_file};
+
+ // utility lambda functions to bind the vector meson and decay particle
+ // types
+ auto momenta_from_tracking =
+ [decay_mass](const std::vector<eic::TrackParametersData>& tracks) {
+ return util::momenta_from_tracking(tracks, decay_mass);
+ };
+
+ //====================================================================
+
+ // Define analysis flow
+ auto d_im =
+ d.Define("p_rec", momenta_from_tracking, {"outputTrackParameters"})
+ .Define("N", "p_rec.size()")
+ .Define("p_sim", util::momenta_from_simulation, {"mcparticles2"})
+ //================================================================
+ .Define("invariant_quantities", util::calc_inv_quant_simu, {"p_sim"})
+ .Define("nu_sim" , util::get_nu_simu, {"invariant_quantities"})
+ .Define("Q2_sim" , util::get_Q2_simu, {"invariant_quantities"})
+ .Define("x_sim" , util::get_x_simu, {"invariant_quantities"});
+ //================================================================
+
+ // Define output histograms
+
+ auto h_nu_sim = d_im.Histo1D(
+ {"h_nu_sim", ";#nu/1000;#", 100, 0., 2.}, "nu_sim");
+ auto h_Q2_sim = d_im.Histo1D(
+ {"h_Q2_sim", ";Q^{2};#", 100, 0., 15.}, "Q2_sim");
+ auto h_x_sim = d_im.Histo1D(
+ {"h_x_sim", ";x;#", 100, 0., 0.1}, "x_sim");
+
+
+ // Plot our histograms.
+ // TODO: to start I'm explicitly plotting the histograms, but want to
+ // factorize out the plotting code moving forward.
+ {
+
+ // Print canvas to output file
+
+ TCanvas c{"canvas2", "canvas2", 1800, 600};
+ c.Divide(3, 1, 0.0001, 0.0001);
+ //pad 1 nu
+ c.cd(1);
+ //gPad->SetLogx(false);
+ //gPad->SetLogy(false);
+ auto& hnu = *h_nu_sim;
+ // histogram style
+ hnu.SetLineColor(plot::kMpBlue);
+ hnu.SetLineWidth(2);
+ // axes
+ hnu.GetXaxis()->CenterTitle();
+ //hnu.GetXaxis()->SetTitle("#times1000");
+ // draw everything
+ hnu.DrawClone("hist");
+ // FIXME hardcoded beam configuration
+ plot::draw_label(10, 100, detector, vm_name, "#nu");
+ TText* tptr21;
+ auto t21 = new TPaveText(.6, .8417, .9, .925, "NB NDC");
+ t21->SetFillColorAlpha(kWhite, 0);
+ t21->SetTextFont(43);
+ t21->SetTextSize(25);
+ tptr21 = t21->AddText("simulated");
+ tptr21->SetTextColor(plot::kMpBlue);
+ //tptr1 = t1->AddText("reconstructed");
+ //tptr1->SetTextColor(plot::kMpOrange);
+ t21->Draw();
+
+ //pad 2 Q2
+ c.cd(2);
+ //gPad->SetLogx(false);
+ //gPad->SetLogy(false);
+ auto& hQ2 = *h_Q2_sim;
+ // histogram style
+ hQ2.SetLineColor(plot::kMpBlue);
+ hQ2.SetLineWidth(2);
+ // axes
+ hQ2.GetXaxis()->CenterTitle();
+ // draw everything
+ hQ2.DrawClone("hist");
+ // FIXME hardcoded beam configuration
+ plot::draw_label(10, 100, detector, vm_name, "Q^{2}");
+ TText* tptr22;
+ auto t22 = new TPaveText(.6, .8417, .9, .925, "NB NDC");
+ t22->SetFillColorAlpha(kWhite, 0);
+ t22->SetTextFont(43);
+ t22->SetTextSize(25);
+ tptr22 = t22->AddText("simulated");
+ tptr22->SetTextColor(plot::kMpBlue);
+ //tptr1 = t1->AddText("reconstructed");
+ //tptr1->SetTextColor(plot::kMpOrange);
+ t22->Draw();
+
+ //pad 1 nu
+ c.cd(3);
+ //gPad->SetLogx(false);
+ //gPad->SetLogy(false);
+ auto& hx = *h_x_sim;
+ // histogram style
+ hx.SetLineColor(plot::kMpBlue);
+ hx.SetLineWidth(2);
+ // axes
+ hx.GetXaxis()->CenterTitle();
+ // draw everything
+ hx.DrawClone("hist");
+ // FIXME hardcoded beam configuration
+ plot::draw_label(10, 100, detector, vm_name, "x");
+ TText* tptr23;
+ auto t23 = new TPaveText(.6, .8417, .9, .925, "NB NDC");
+ t23->SetFillColorAlpha(kWhite, 0);
+ t23->SetTextFont(43);
+ t23->SetTextSize(25);
+ tptr23 = t23->AddText("simulated");
+ tptr23->SetTextColor(plot::kMpBlue);
+ //tptr1 = t1->AddText("reconstructed");
+ //tptr1->SetTextColor(plot::kMpOrange);
+ t23->Draw();
+
+ c.Print(fmt::format("{}InvariantQuantities.png", output_prefix).c_str());
+ }
+
+ // TODO we're not actually doing an IM fit yet, so for now just return an
+ // error for the test result
+ vm_mass_resolution_test.error(-1);
+
+ // write out our test data
+ eic::util::write_test(vm_mass_resolution_test,
+ fmt::format("{}vm_invar.json", output_prefix));
+
+ // That's all!
+ return 0;
+}
diff --git a/dvmp/analysis/vm_mass.cxx b/dvmp/analysis/vm_mass.cxx
index 70de8261..fb736909 100644
--- a/dvmp/analysis/vm_mass.cxx
+++ b/dvmp/analysis/vm_mass.cxx
@@ -19,6 +19,7 @@
// file prefix), and labeled with our detector name.
// TODO: I think it would be better to pass small json configuration file to
// the test, instead of this ever-expanding list of function arguments.
+// FIXME: MC does not trace back into particle history. Need to fix that
int vm_mass(const std::string& config_name) {
// read our configuration
std::ifstream config_file{config_name};
@@ -77,7 +78,9 @@ int vm_mass(const std::string& config_name) {
[vm_mass](const std::vector<ROOT::Math::PxPyPzMVector>& parts) {
return util::find_decay_pair(parts, vm_mass);
};
-
+
+
+ //util::PrintGeant4(mcparticles2);
// Define analysis flow
auto d_im =
d.Define("p_rec", momenta_from_tracking, {"outputTrackParameters"})
@@ -93,6 +96,7 @@ int vm_mass(const std::string& config_name) {
.Define("phi_sim" , util::get_phi, {"decay_pair_sim"})
.Define("rapidity_rec" , util::get_y, {"decay_pair_rec"})
.Define("rapidity_sim" , util::get_y, {"decay_pair_sim"});
+
// Define output histograms
auto h_im_rec = d_im.Histo1D(
@@ -116,6 +120,7 @@ int vm_mass(const std::string& config_name) {
{"h_y_sim", ";y_{ll'};#", 1000, -5., 5.}, "rapidity_sim");
+
// Plot our histograms.
// TODO: to start I'm explicitly plotting the histograms, but want to
// factorize out the plotting code moving forward.
@@ -241,9 +246,9 @@ int vm_mass(const std::string& config_name) {
tptr4 = t4->AddText("reconstructed");
tptr4->SetTextColor(plot::kMpOrange);
t4->Draw();
-
- // Print canvas to output file
+
c.Print(fmt::format("{}vm_mass_pt_phi_rapidity.png", output_prefix).c_str());
+
}
// TODO we're not actually doing an IM fit yet, so for now just return an
diff --git a/dvmp/dvmp.sh b/dvmp/dvmp.sh
index acad1d1c..6b1735ec 100755
--- a/dvmp/dvmp.sh
+++ b/dvmp/dvmp.sh
@@ -56,6 +56,7 @@ GEN_FILE=${INPUT_PATH}/gen-${CONFIG}_${DECAY}.hepmc
SIM_FILE=${TMP_PATH}/sim-${CONFIG}_${DECAY}.root
SIM_LOG=${TMP_PATH}/sim-${CONFIG}_${DECAY}.log
+
REC_FILE=${TMP_PATH}/rec-${CONFIG}_${DECAY}.root
REC_LOG=${TMP_PATH}/sim-${CONFIG}_${DECAY}.log
@@ -65,7 +66,7 @@ PLOT_TAG=${CONFIG}_${DECAY}
## Step 2: Run the simulation
echo "Running Geant4 simulation"
npsim --runType batch \
- --part.minimalKineticEnergy 1000*GeV \
+ --part.minimalKineticEnergy 100*GeV \
-v WARNING \
--numberOfEvents ${JUGGLER_N_EVENTS} \
--compactFile ${DETECTOR_PATH}/${JUGGLER_DETECTOR}.xml \
@@ -127,7 +128,7 @@ EOF
## run the analysis script with this configuration
root -b -q "dvmp/analysis/vm_mass.cxx(\"${CONFIG}\")"
-
+root -b -q "dvmp/analysis/vm_invar.cxx(\"${CONFIG}\")"
if [ "$?" -ne "0" ] ; then
echo "ERROR running root script"
exit 1
@@ -146,6 +147,7 @@ fi
## Always move over log files to the results path
mv ${SIM_LOG} ${REC_LOG} ${RESULTS_PATH}
+
## cleanup output files
rm -f ${REC_FILE} ${SIM_FILE}
--
GitLab