Skip to content
Snippets Groups Projects

Add y, Q2, x, t resolution tests

Open Add y, Q2, x, t resolution tests
ziyue_work_branch into master
Open Ziyue Zhang requested to merge ziyue_work_branch into master
Compare
and
4 files
+ 160
− 79
Compare changes
  • Side-by-side
  • Inline
Files
4
benchmarks/dvmp/analysis/vm_invar.cxx
+ 147
− 71
@@ -43,16 +43,46 @@ int vm_invar(const std::string& config_name)
@@ -43,16 +43,46 @@ int vm_invar(const std::string& config_name)
fmt::print(" - output prefix: {}\n", output_prefix);
fmt::print(" - output prefix: {}\n", output_prefix);
// create our test definition
// create our test definition
// test_tag
std::vector<eic::util::Test> Tests;
 
eic::util::Test y_resolution_test{
 
{{"name", fmt::format("{}_y_resolution", test_tag)},
 
{"title",
 
fmt::format("y Resolution for {} -> {} events with {}", vm_name, decay_name, detector)},
 
{"description", "y resolution: relative difference with Gaussian fit"},
 
{"quantity", "resolution"},
 
{"target", ".4"}}};
 
Tests.push_back(y_resolution_test);
 
eic::util::Test Q2_resolution_test{
eic::util::Test Q2_resolution_test{
{{"name", fmt::format("{}_Q2_resolution", test_tag)},
{{"name", fmt::format("{}_Q2_resolution", test_tag)},
{"title",
{"title",
fmt::format("Q^2 Resolution for {} -> {} events with {}", vm_name, decay_name, detector)},
fmt::format("Q^2 Resolution for {} -> {} events with {}", vm_name, decay_name, detector)},
{"description", "Invariant Mass Resolution calculated from raw "
{"description", "Q^2 resolution: relative difference with Gaussian fit"},
"tracking data using a Gaussian fit."},
{"quantity", "resolution"},
 
{"target", ".09"}}};
 
Tests.push_back(Q2_resolution_test);
 
 
eic::util::Test x_resolution_test{
 
{{"name", fmt::format("{}_x_resolution", test_tag)},
 
{"title",
 
fmt::format("x Resolution for {} -> {} events with {}", vm_name, decay_name, detector)},
 
{"description", "x resolution: relative difference with Gaussian fit"},
 
{"quantity", "resolution"},
 
{"target", ".35"}}};
 
Tests.push_back(x_resolution_test);
 
 
eic::util::Test t_resolution_test{
 
{{"name", fmt::format("{}_t_resolution", test_tag)},
 
{"title",
 
fmt::format("t Resolution for {} -> {} events with {}", vm_name, decay_name, detector)},
 
{"description", "t resolution: relative difference with Gaussian fit"},
{"quantity", "resolution"},
{"quantity", "resolution"},
{"target", ".1"}}};
{"target", ".07"}}};
 
Tests.push_back(t_resolution_test);
 
double width_target[4] = {0.4, 0.09, 0.35, 0.07};
 
TH1::SetDefaultSumw2();
 
// Run this in multi-threaded mode if desired
// Run this in multi-threaded mode if desired
ROOT::EnableImplicitMT(kNumThreads);
ROOT::EnableImplicitMT(kNumThreads);
@@ -92,48 +122,103 @@ int vm_invar(const std::string& config_name)
@@ -92,48 +122,103 @@ int vm_invar(const std::string& config_name)
.Define("y_sim", util::get_y, {"invariant_quantities_sim"})
.Define("y_sim", util::get_y, {"invariant_quantities_sim"})
.Define("Q2_sim", util::get_Q2, {"invariant_quantities_sim"})
.Define("Q2_sim", util::get_Q2, {"invariant_quantities_sim"})
.Define("x_sim", util::get_x, {"invariant_quantities_sim"})
.Define("x_sim", util::get_x, {"invariant_quantities_sim"})
.Define("t_sim", util::get_t, {"invariant_quantities_sim"});
.Define("t_sim", util::get_t, {"invariant_quantities_sim"})
 
.Define("y_diff", "(y_rec - y_sim)/y_sim")
 
.Define("Q2_diff", "(Q2_rec - Q2_sim)/Q2_sim")
 
.Define("x_diff", "(x_rec - x_sim)/x_sim")
 
.Define("t_diff", "(t_rec - t_sim)/t_sim");
//================================================================
//================================================================
 
//Factorizeation
 
double func_range[4] = {1.5, 0.3, 1., 2.};
 
double hist_range_l[4] = {0., 0., 0., -1.};
 
double hist_range_h[4] = {1., 15., 0.1, 0.};
 
 
/*std::string VarName[4] = {"y", "Q2", "x", "t"};
 
std::string histName[4];
 
std::string histTitle[4];
 
std::string RawHist[4];
 
for(int i = 0 ; i < 4 ; i++){
 
histName[i] = "h_" + VarName[i] + "_sim_test";
 
if(i!=1){
 
histTitle[i] = ";" + VarName[i] + ";#";
 
}else{
 
histTitle[i] = ";Q^{2};#";
 
}
 
RawHist[i] = VarName[i] + "_sim";
 
}
 
 
TH1D* h_sim[4];
 
{
 
int i = 0;
 
cout<<"================"<<histName[i]<<"================"<<endl;
 
//auto h_tmp = d_im.Histo1D({histName[i], ";y;#", 50, hist_range_l[i], hist_range_h[i]}, "y_sim"); //using string variable
 
//auto h_tmp = d_im.Histo1D({"h_y_sim_test", ";y;#", 50, hist_range_l[i], hist_range_h[i]}, "y_sim"); //directly quote the string
 
}*/
 
//==================================================================
// Define output histograms
// Define output histograms
//auto h_nu_sim = d_im.Histo1D({"h_nu_sim", ";#nu/1000;#", 100, 0., 2.}, "nu_sim");
//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};#", 50, 0., 15.}, "Q2_sim");
auto h_Q2_sim = d_im.Histo1D({"h_Q2_sim", ";Q^{2};#", 50, 0., 15.}, "Q2_sim");
auto h_x_sim = d_im.Histo1D({"h_x_sim", ";x;#", 50, 0., 0.1}, "x_sim");
auto h_x_sim = d_im.Histo1D({"h_x_sim", ";x;#", 50, 0., 0.1}, "x_sim");
auto h_y_sim = d_im.Histo1D({"h_y_sim", ";y;#", 50, 0., 1.}, "y_sim");
auto h_y_sim = d_im.Histo1D({"h_y_sim", ";y;#", 50, 0., 1.}, "y_sim");
auto h_t_sim = d_im.Histo1D({"h_t_sim", ";t;#", 50, -1., 0.}, "t_sim");
auto h_t_sim = d_im.Histo1D({"h_t_sim", ";t;#", 50, -1., 0.}, "t_sim");
//auto h_nu_rec = d_im.Histo1D({"h_nu_rec", ";#nu/1000;#", 100, 0., 2.}, "nu_rec");
//auto h_nu_rec = d_im.Histo1D({"h_nu_rec", ";#nu/1000;#", 100, 0., 2.}, "nu_rec");
auto h_Q2_rec = d_im.Histo1D({"h_Q2_rec", ";Q^{2};#", 50, 0., 15.}, "Q2_rec");
auto h_Q2_rec = d_im.Histo1D({"h_Q2_rec", ";Q^{2};#", 50, 0., 15.}, "Q2_rec");
auto h_x_rec = d_im.Histo1D({"h_x_rec", ";x;#", 50, 0., 0.1}, "x_rec");
auto h_x_rec = d_im.Histo1D({"h_x_rec", ";x;#", 50, 0., 0.1}, "x_rec");
auto h_y_rec = d_im.Histo1D({"h_y_rec", ";y;#", 50, 0., 1.}, "y_rec");
auto h_y_rec = d_im.Histo1D({"h_y_rec", ";y;#", 50, 0., 1.}, "y_rec");
auto h_t_rec = d_im.Histo1D({"h_t_rec", ";t;#", 50, -1., 0.}, "t_rec");
auto h_t_rec = d_im.Histo1D({"h_t_rec", ";t;#", 50, -1., 0.}, "t_rec");
 
auto h_y_diff = d_im.Histo1D({"h_y_diff", ";#Deltay/y;#", 50, -1.5, 1.5}, "y_diff");
 
auto h_Q2_diff = d_im.Histo1D({"h_Q2_diff", ";#DeltaQ^{2}/Q^{2};#", 50, -0.3, 0.3}, "Q2_diff");
 
auto h_x_diff = d_im.Histo1D({"h_x_diff", ";#Deltax/x;#", 50, -1., 1.}, "x_diff");
 
auto h_t_diff = d_im.Histo1D({"h_t_diff", ";#Deltat/t;#", 50, -0.5, 0.5}, "t_diff");
 
 
double nEvents = h_y_diff->Integral(0, -1);
 
 
TH1D* histtest = (TH1D*)d_im.Histo1D({"h_Q2_sim_test", ";Q^{2};#", 50, 0., 15.}, "Q2_sim");
 
TCanvas c{"ctest", "ctest", 1200, 900};
 
histtest->Draw("hist e");
 
c.Print(fmt::format("{}test.png", output_prefix).c_str());
 
// 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
// factorize out the plotting code moving forward.
// factorize out the plotting code moving forward.
{
TFitResultPtr myFitPtr[4];
TF1* myf[4];
 
for(int i = 0 ; i < 4 ; i++){
 
myf[i] = new TF1(Form("myf_%d", i), "[2]*TMath::Gaus(x, [0], [1], 0)", -func_range[i], func_range[i]);
 
myf[i]->SetParameters(0., 0.25, nEvents/10.);
 
myf[i]->SetParLimits(0, -0.5, 0.5);
 
myf[i]->SetParLimits(1, 0., 1.0);
 
/*if(i==3){
 
myf[i]->SetParameter(1, 0.1);
 
//myf[i]->SetParameter(2, nEvents);
 
}*/
 
myf[i]->SetParLimits(2, 0., nEvents*10.);
 
myf[i]->SetNpx(1000);
 
myf[i]->SetLineColor(2);
 
myf[i]->SetLineStyle(7);
 
}
 
 
// Print canvas to output file
// Print canvas to output file
TCanvas c{"canvas2", "canvas2", 1200, 900};
TCanvas c{"canvas2", "canvas2", 1200, 900};
c.Divide(2, 2, 0.0001, 0.0001);
c.Divide(2, 2, 0.0001, 0.0001);
// pad 1 nu
//============================================================================
 
//pad 1 nu_diff
c.cd(1);
c.cd(1);
auto& hy_rec = *h_y_rec;
auto& hy_diff = *h_y_diff;
auto& hy_sim = *h_y_sim;
// histogram style
// histogram style
hy_rec.SetLineColor(plot::kMpOrange);
hy_diff.SetLineColor(plot::kMpOrange);
hy_rec.SetLineWidth(1);
hy_diff.SetLineWidth(1);
hy_sim.SetLineColor(plot::kMpBlue);
hy_sim.SetLineWidth(2);
// axes
// axes
hy_sim.GetXaxis()->CenterTitle();
hy_diff.GetXaxis()->CenterTitle();
// draw everything
// draw everything
hy_sim.DrawClone("hist");
hy_diff.DrawClone("hist");
hy_rec.DrawClone("hist same");
myFitPtr[0] = hy_diff.Fit(myf[0], "S 0", "", -func_range[0], func_range[0]);
 
myf[0]->Draw("same");
// FIXME hardcoded beam configuration
// FIXME hardcoded beam configuration
plot::draw_label(10, 100, detector);
plot::draw_label(10, 100, detector);
TText* tptr1;
TText* tptr1;
@@ -141,26 +226,23 @@ int vm_invar(const std::string& config_name)
@@ -141,26 +226,23 @@ int vm_invar(const std::string& config_name)
t1->SetFillColorAlpha(kWhite, 0);
t1->SetFillColorAlpha(kWhite, 0);
t1->SetTextFont(43);
t1->SetTextFont(43);
t1->SetTextSize(25);
t1->SetTextSize(25);
tptr1 = t1->AddText("simulated");
tptr1 = t1->AddText("#Deltay/y");
tptr1->SetTextColor(plot::kMpBlue);
tptr1 = t1->AddText("rec(PlaceHolder)");
tptr1->SetTextColor(plot::kMpOrange);
tptr1->SetTextColor(plot::kMpOrange);
t1->Draw();
t1->Draw();
// pad 2 Q2
 
// pad 2 Q2_diff
c.cd(2);
c.cd(2);
auto& hQ2_rec = *h_Q2_rec;
auto& hQ2_diff = *h_Q2_diff;
auto& hQ2_sim = *h_Q2_sim;
// histogram style
// histogram style
hQ2_rec.SetLineColor(plot::kMpOrange);
hQ2_diff.SetLineColor(plot::kMpOrange);
hQ2_rec.SetLineWidth(1);
hQ2_diff.SetLineWidth(1);
hQ2_sim.SetLineColor(plot::kMpBlue);
hQ2_sim.SetLineWidth(2);
// axes
// axes
hQ2_sim.GetXaxis()->CenterTitle();
hQ2_diff.GetXaxis()->CenterTitle();
// draw everything
// draw everything
hQ2_sim.DrawClone("hist");
hQ2_diff.DrawClone("hist");
hQ2_rec.DrawClone("hist same");
myFitPtr[1] = hQ2_diff.Fit(myf[1], "S 0", "", -func_range[1], func_range[1]);
 
myf[1]->Draw("same");
// FIXME hardcoded beam configuration
// FIXME hardcoded beam configuration
plot::draw_label(10, 100, detector);
plot::draw_label(10, 100, detector);
TText* tptr2;
TText* tptr2;
@@ -168,26 +250,22 @@ int vm_invar(const std::string& config_name)
@@ -168,26 +250,22 @@ int vm_invar(const std::string& config_name)
t2->SetFillColorAlpha(kWhite, 0);
t2->SetFillColorAlpha(kWhite, 0);
t2->SetTextFont(43);
t2->SetTextFont(43);
t2->SetTextSize(25);
t2->SetTextSize(25);
tptr2 = t2->AddText("simulated");
tptr2 = t2->AddText("#DeltaQ^{2}/Q^{2}");
tptr2->SetTextColor(plot::kMpBlue);
tptr2 = t2->AddText("rec(PlaceHolder)");
tptr2->SetTextColor(plot::kMpOrange);
tptr2->SetTextColor(plot::kMpOrange);
t2->Draw();
t2->Draw();
// pad 3 x
// pad 3 x_diff
c.cd(3);
c.cd(3);
auto& hx_rec = *h_x_rec;
auto& hx_diff = *h_x_diff;
auto& hx_sim = *h_x_sim;
// histogram style
// histogram style
hx_rec.SetLineColor(plot::kMpOrange);
hx_diff.SetLineColor(plot::kMpOrange);
hx_rec.SetLineWidth(1);
hx_diff.SetLineWidth(1);
hx_sim.SetLineColor(plot::kMpBlue);
hx_sim.SetLineWidth(2);
// axes
// axes
hx_sim.GetXaxis()->CenterTitle();
hx_diff.GetXaxis()->CenterTitle();
// draw everything
// draw everything
hx_sim.DrawClone("hist");
hx_diff.DrawClone("hist");
hx_rec.DrawClone("hist same");
myFitPtr[2] = hx_diff.Fit(myf[2], "S 0", "", -func_range[2], func_range[2]);
 
myf[2]->Draw("same");
// FIXME hardcoded beam configuration
// FIXME hardcoded beam configuration
plot::draw_label(10, 100, detector);
plot::draw_label(10, 100, detector);
TText* tptr3;
TText* tptr3;
@@ -195,26 +273,22 @@ int vm_invar(const std::string& config_name)
@@ -195,26 +273,22 @@ int vm_invar(const std::string& config_name)
t3->SetFillColorAlpha(kWhite, 0);
t3->SetFillColorAlpha(kWhite, 0);
t3->SetTextFont(43);
t3->SetTextFont(43);
t3->SetTextSize(25);
t3->SetTextSize(25);
tptr3 = t3->AddText("simulated");
tptr3 = t3->AddText("#Deltax/x");
tptr3->SetTextColor(plot::kMpBlue);
tptr3 = t3->AddText("rec(PlaceHolder)");
tptr3->SetTextColor(plot::kMpOrange);
tptr3->SetTextColor(plot::kMpOrange);
t3->Draw();
t3->Draw();
// pad 4 t
// pad 4 t_diff
c.cd(4);
c.cd(4);
auto& ht_rec = *h_t_rec;
auto& ht_diff = *h_t_diff;
auto& ht_sim = *h_t_sim;
// histogram style
// histogram style
ht_rec.SetLineColor(plot::kMpOrange);
ht_diff.SetLineColor(plot::kMpOrange);
ht_rec.SetLineWidth(1);
ht_diff.SetLineWidth(1);
ht_sim.SetLineColor(plot::kMpBlue);
ht_sim.SetLineWidth(2);
// axes
// axes
ht_sim.GetXaxis()->CenterTitle();
ht_diff.GetXaxis()->CenterTitle();
// draw everything
// draw everything
ht_sim.DrawClone("hist");
ht_diff.DrawClone("hist");
ht_rec.DrawClone("hist same");
myFitPtr[3] = ht_diff.Fit(myf[3], "S 0", "", -func_range[3], func_range[3]);
 
myf[3]->Draw("same");
// FIXME hardcoded beam configuration
// FIXME hardcoded beam configuration
plot::draw_label(10, 100, detector);
plot::draw_label(10, 100, detector);
TText* tptr4;
TText* tptr4;
@@ -222,23 +296,25 @@ int vm_invar(const std::string& config_name)
@@ -222,23 +296,25 @@ int vm_invar(const std::string& config_name)
t4->SetFillColorAlpha(kWhite, 0);
t4->SetFillColorAlpha(kWhite, 0);
t4->SetTextFont(43);
t4->SetTextFont(43);
t4->SetTextSize(25);
t4->SetTextSize(25);
tptr4 = t4->AddText("simulated");
tptr4 = t4->AddText("#Deltat/t");
tptr4->SetTextColor(plot::kMpBlue);
tptr4 = t4->AddText("rec(PlaceHolder)");
tptr4->SetTextColor(plot::kMpOrange);
tptr4->SetTextColor(plot::kMpOrange);
t4->Draw();
t4->Draw();
//============================================================================
c.Print(fmt::format("{}InvariantQuantities.png", output_prefix).c_str());
c.Print(fmt::format("{}InvariantQuantities.png", output_prefix).c_str());
}
// TODO we're not actually getting the resolutions yet
// error for the test result
Q2_resolution_test.error(-1);
 
for(int i = 0 ; i < 4 ; i++){
 
double width = myf[i]->GetParameter(1);
 
if(myFitPtr[i]->Status()!=0){
 
Tests[i].error(-1);
 
}else if(width > width_target[i]){
 
Tests[i].fail(width);
 
}else{
 
Tests[i].pass(width);
 
}
 
}
 
// write out our test data
// write out our test data
eic::util::write_test(Q2_resolution_test, fmt::format("{}invar.json", output_prefix));
eic::util::write_test(Tests, fmt::format("{}invar.json", output_prefix));
// That's all!
// That's all!
return 0;
return 0;