Skip to content
Snippets Groups Projects

Add y, Q2, x, t resolution tests

Open Ziyue Zhang requested to merge ziyue_work_branch into master
Compare and
3 files
+ 193
162
Compare changes
  • Side-by-side
  • Inline
Files
3
@@ -41,21 +41,27 @@ int vm_invar(const std::string& config_name)
fmt::print(" - Detector package: {}\n", detector);
fmt::print(" - input file: {}\n", rec_file);
fmt::print(" - output prefix: {}\n", output_prefix);
// create our test definition
// test_tag
eic::util::Test Q2_resolution_test{
{{"name", fmt::format("{}_Q2_resolution", test_tag)},
std::string VarName[4] = {"y", "Q2", "x", "t"};
double width_target[4] = {.42, .093, .34, .071};
//============================== test definition ==============================
std::vector<eic::util::Test> Tests;
for(int i = 0 ; i < 4 ; i++){
eic::util::Test resolution_test_tmp{
{{"name", fmt::format("{}_{}_resolution", test_tag, VarName[i])},
{"title",
fmt::format("Q^2 Resolution for {} -> {} events with {}", vm_name, decay_name, detector)},
{"description", "Invariant Mass Resolution calculated from raw "
"tracking data using a Gaussian fit."},
fmt::format("{} Resolution for {} -> {} events with {}", VarName[i], vm_name, decay_name, detector)},
{"description", fmt::format("{} resolution: relative difference with Gaussian fit", VarName[i])},
{"quantity", "resolution"},
{"target", ".1"}}};
{"target", fmt::format("{}", width_target[i])}}};
Tests.push_back(resolution_test_tmp);
}
//============================== general settings ==============================
// Run this in multi-threaded mode if desired
ROOT::EnableImplicitMT(kNumThreads);
TH1::SetDefaultSumw2();
// 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);
@@ -67,19 +73,16 @@ int vm_invar(const std::string& config_name)
// 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
//============================== redef functions ==============================
auto momenta_sort_sim = [vm_name, decay_name](const std::vector<dd4pod::Geant4ParticleData>& parts){
return util::momenta_sort_sim(parts, vm_name, decay_name);
};
auto momenta_sort_rec = [vm_name, decay_name](const std::vector<eic::ReconstructedParticleData>& parts){
return util::momenta_sort_rec(parts, vm_name, decay_name);
};
//====================================================================
// Define analysis flow
//============================== Define analysis flow ==============================
auto d_im = d.Define("p_rec_sorted", momenta_sort_rec, {"DummyReconstructedParticles"})
.Define("p_sim_sorted", momenta_sort_sim, {"mcparticles2"})
.Define("N", "p_rec_sorted.size()")
@@ -92,154 +95,177 @@ int vm_invar(const std::string& config_name)
.Define("y_sim", util::get_y, {"invariant_quantities_sim"})
.Define("Q2_sim", util::get_Q2, {"invariant_quantities_sim"})
.Define("x_sim", util::get_x, {"invariant_quantities_sim"})
.Define("t_sim", util::get_t, {"invariant_quantities_sim"});
//================================================================
// 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};#", 50, 0., 15.}, "Q2_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_t_sim = d_im.Histo1D({"h_t_sim", ";t;#", 50, -1., 0.}, "t_sim");
.Define("t_sim", util::get_t, {"invariant_quantities_sim"})
.Define("y_dif", "y_rec - y_sim")
.Define("Q2_dif", "Q2_rec - Q2_sim")
.Define("x_dif", "x_rec - x_sim")
.Define("t_dif", "t_rec - t_sim")
.Define("y_rdf", "(y_rec - y_sim)/y_sim")
.Define("Q2_rdf", "(Q2_rec - Q2_sim)/Q2_sim")
.Define("x_rdf", "(x_rec - x_sim)/x_sim")
.Define("t_rdf", "(t_rec - t_sim)/t_sim");
//==============================hist def==============================
//ranges
double range_l[4][4] = {{0., 0., -0.3, -1.5}, { 0., 0., -1., -0.3}, {0.0, 0.0, -0.01, -1.}, {-1., -1., -0.1, -0.5}};
double range_h[4][4] = {{1., 1., 0.3, 1.5}, {15., 15., 1., 0.3}, {0.1, 0.1, 0.01, 1.}, { 0., 0., 0.1, 0.5}};
//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_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_t_rec = d_im.Histo1D({"h_t_rec", ";t;#", 50, -1., 0.}, "t_rec");
//strings used
std::string VarCate[5] = {"sim", "rec", "dif", "rdf", "svr"};
std::string histName[4][5];
std::string histTitles[4][5];
std::string RawhistName[4][4];
// 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", 1200, 900};
c.Divide(2, 2, 0.0001, 0.0001);
// pad 1 nu
c.cd(1);
auto& hy_rec = *h_y_rec;
auto& hy_sim = *h_y_sim;
// histogram style
hy_rec.SetLineColor(plot::kMpOrange);
hy_rec.SetLineWidth(1);
hy_sim.SetLineColor(plot::kMpBlue);
hy_sim.SetLineWidth(2);
// axes
hy_sim.GetXaxis()->CenterTitle();
// draw everything
hy_sim.DrawClone("hist");
hy_rec.DrawClone("hist same");
// FIXME hardcoded beam configuration
plot::draw_label(10, 100, detector);
TText* tptr1;
auto t1 = new TPaveText(.6, .8417, .9, .925, "NB NDC");
t1->SetFillColorAlpha(kWhite, 0);
t1->SetTextFont(43);
t1->SetTextSize(25);
tptr1 = t1->AddText("simulated");
tptr1->SetTextColor(plot::kMpBlue);
tptr1 = t1->AddText("rec(PlaceHolder)");
tptr1->SetTextColor(plot::kMpOrange);
t1->Draw();
// pad 2 Q2
c.cd(2);
auto& hQ2_rec = *h_Q2_rec;
auto& hQ2_sim = *h_Q2_sim;
// histogram style
hQ2_rec.SetLineColor(plot::kMpOrange);
hQ2_rec.SetLineWidth(1);
hQ2_sim.SetLineColor(plot::kMpBlue);
hQ2_sim.SetLineWidth(2);
// axes
hQ2_sim.GetXaxis()->CenterTitle();
// draw everything
hQ2_sim.DrawClone("hist");
hQ2_rec.DrawClone("hist same");
// FIXME hardcoded beam configuration
plot::draw_label(10, 100, detector);
TText* tptr2;
auto t2 = new TPaveText(.6, .8417, .9, .925, "NB NDC");
t2->SetFillColorAlpha(kWhite, 0);
t2->SetTextFont(43);
t2->SetTextSize(25);
tptr2 = t2->AddText("simulated");
tptr2->SetTextColor(plot::kMpBlue);
tptr2 = t2->AddText("rec(PlaceHolder)");
tptr2->SetTextColor(plot::kMpOrange);
t2->Draw();
TH1D* h_Var1D[4][4];
TH2D* h_Var2D[4];
for(int i = 0 ; i < 4 ; i++){
for(int j = 0 ; j < 4 ; j++){
//construct histName
histName[i][j] = "h_" + VarName[i] + "_" + VarCate[j];
//construct histTitles
histTitles[i][j] = ";";
if(j > 1) histTitles[i][j] = histTitles[i][j] + "#Delta";
if(i==1){
histTitles[i][j] = histTitles[i][j] + "Q^{2}";
}else{
histTitles[i][j] = histTitles[i][j] + VarName[i];
}
if(j==3){
histTitles[i][j] = histTitles[i][j] + "/";
if(i==1){
histTitles[i][j] = histTitles[i][j] + "Q^{2}";
}else{
histTitles[i][j] = histTitles[i][j] + VarName[i];
}
}
histTitles[i][j] = histTitles[i][j] + ";#";
//construct RawhistName
RawhistName[i][j] = VarName[i] + "_" + VarCate[j];
//get histograms
auto h_tmp = d_im.Histo1D({fmt::format("{}_tmp", histName[i][j]).c_str(), histTitles[i][j].c_str(), 50, range_l[i][j], range_h[i][j]}, RawhistName[i][j].c_str());
auto& htmp = *h_tmp;
h_Var1D[i][j] = (TH1D*)htmp.Clone(histName[i][j].c_str());
}
//2d histogram
histName[i][4] = "h_" + VarName[i] + "_" + VarCate[4];
if(i==1){
histTitles[i][4] = ";Q^{2}_{sim};Q^{2}_{rec};#";
}else{
histTitles[i][4] = ";" + VarName[i] + "_{sim};" + VarName[i] + "_{rec};#";
}
auto h_tmp = d_im.Histo2D({fmt::format("{}_tmp", histName[i][4]).c_str(), histTitles[i][4].c_str(), 50, range_l[i][0], range_h[i][0], 50, range_l[i][0], range_h[i][0]}, RawhistName[i][0].c_str(), RawhistName[i][1].c_str());
auto& htmp = *h_tmp;
h_Var2D[i] = (TH2D*)htmp.Clone(histName[i][4].c_str());
}
double nEvents = h_Var1D[0][0]->Integral(0, -1);
//==============================fit and plot==============================
TFitResultPtr myFitPtr[4][2];
TF1* myf[4][2];
TText* tptr[4][4];
TPaveText* t[4][4];
for(int i = 0 ; i < 4 ; i++){
TCanvas* ctmp = new TCanvas("ctmp", "ctmp", 1200, 900);
ctmp->Divide(2, 2, 0.001, 0.001);
//for pad 1: sim rec overlay
h_Var1D[i][0]->SetLineColor(plot::kMpBlue);
h_Var1D[i][0]->SetLineWidth(2);
h_Var1D[i][0]->GetXaxis()->CenterTitle();
h_Var1D[i][1]->SetLineColor(plot::kMpOrange);
h_Var1D[i][1]->SetLineWidth(1);
// pad 3 x
c.cd(3);
auto& hx_rec = *h_x_rec;
auto& hx_sim = *h_x_sim;
// histogram style
hx_rec.SetLineColor(plot::kMpOrange);
hx_rec.SetLineWidth(1);
hx_sim.SetLineColor(plot::kMpBlue);
hx_sim.SetLineWidth(2);
// axes
hx_sim.GetXaxis()->CenterTitle();
// draw everything
hx_sim.DrawClone("hist");
hx_rec.DrawClone("hist same");
// FIXME hardcoded beam configuration
plot::draw_label(10, 100, detector);
TText* tptr3;
auto t3 = new TPaveText(.6, .8417, .9, .925, "NB NDC");
t3->SetFillColorAlpha(kWhite, 0);
t3->SetTextFont(43);
t3->SetTextSize(25);
tptr3 = t3->AddText("simulated");
tptr3->SetTextColor(plot::kMpBlue);
tptr3 = t3->AddText("rec(PlaceHolder)");
tptr3->SetTextColor(plot::kMpOrange);
t3->Draw();
//for pad 2: sim~rec 2d
h_Var2D[i]->GetXaxis()->CenterTitle();
h_Var2D[i]->GetYaxis()->CenterTitle();
// pad 4 t
c.cd(4);
auto& ht_rec = *h_t_rec;
auto& ht_sim = *h_t_sim;
// histogram style
ht_rec.SetLineColor(plot::kMpOrange);
ht_rec.SetLineWidth(1);
ht_sim.SetLineColor(plot::kMpBlue);
ht_sim.SetLineWidth(2);
// axes
ht_sim.GetXaxis()->CenterTitle();
// draw everything
ht_sim.DrawClone("hist");
ht_rec.DrawClone("hist same");
// FIXME hardcoded beam configuration
plot::draw_label(10, 100, detector);
TText* tptr4;
auto t4 = new TPaveText(.6, .8417, .9, .925, "NB NDC");
t4->SetFillColorAlpha(kWhite, 0);
t4->SetTextFont(43);
t4->SetTextSize(25);
tptr4 = t4->AddText("simulated");
tptr4->SetTextColor(plot::kMpBlue);
tptr4 = t4->AddText("rec(PlaceHolder)");
tptr4->SetTextColor(plot::kMpOrange);
t4->Draw();
//for pad 3: rec - sim
h_Var1D[i][2]->SetLineColor(plot::kMpGrey);
h_Var1D[i][2]->SetLineWidth(1);
h_Var1D[i][2]->GetXaxis()->CenterTitle();
c.Print(fmt::format("{}InvariantQuantities.png", output_prefix).c_str());
//for pad 4: (rec - sim)/sim
h_Var1D[i][3]->SetLineColor(plot::kMpGrey);
h_Var1D[i][3]->SetLineWidth(1);
h_Var1D[i][3]->GetXaxis()->CenterTitle();
//initialize myf
for(int j = 0 ; j < 2 ; j++){
myf[i][j] = new TF1(Form("myf_%d_%d", i, j), "[2]*TMath::Gaus(x, [0], [1], 0)", range_l[i][j+2], range_h[i][j+2]);
myf[i][j]->SetParameters(0., (range_h[i][j+2]-range_l[i][j+2])/4., nEvents/20.);
myf[i][j]->SetParLimits(0, -0.5, 0.5);
myf[i][j]->SetParLimits(1, 0., 1.0);
myf[i][j]->SetParLimits(2, 0., nEvents*1000.);
myf[i][j]->SetNpx(1000);
myf[i][j]->SetLineColor(plot::kMpRed);
myf[i][j]->SetLineStyle(7);
}
//fit and plot
for(int j = 0 ; j < 4 ; j++){
ctmp->cd(j+1);
t[i][j] = new TPaveText(.6, .8417, .9, .925, "NB NDC");
t[i][j]->SetFillColorAlpha(kWhite, 0);
t[i][j]->SetTextFont(43);
t[i][j]->SetTextSize(25);
switch(j){
case 0://sim rec overlay
h_Var1D[i][0]->Draw("hist");
h_Var1D[i][1]->Draw("hist same");
tptr[i][j] = t[i][j]->AddText("simulation");
tptr[i][j]->SetTextColor(plot::kMpBlue);
tptr[i][j] = t[i][j]->AddText("reconstructed");
tptr[i][j]->SetTextColor(plot::kMpOrange);
plot::draw_label(10, 100, detector);
break;
case 1:
gPad->SetLogz();
h_Var2D[i]->Draw("colz");
break;//2d
case 2://dx
h_Var1D[i][2]->Draw("hist");
myFitPtr[i][0] = h_Var1D[i][2]->Fit(myf[i][0], "S 0", "", range_l[i][j], range_h[i][j]);
myf[i][0]->Draw("same");
if(i==1){
tptr[i][j] = t[i][j]->AddText("#DeltaQ^{2}");
}else{
tptr[i][j] = t[i][j]->AddText(fmt::format("#Delta{}", VarName[i]).c_str());
}
tptr[i][j]->SetTextColor(plot::kMpGrey);
break;
case 3://dx/x
h_Var1D[i][3]->Draw("hist");
myFitPtr[i][1] = h_Var1D[i][3]->Fit(myf[i][1], "S 0", "", range_l[i][j], range_h[i][j]);
myf[i][1]->Draw("same");
if(i==1){
tptr[i][j] = t[i][j]->AddText("#DeltaQ^{2}/Q^{2}");
}else{
tptr[i][j] = t[i][j]->AddText(fmt::format("#Delta{}/{}", VarName[i], VarName[i]).c_str());
}
tptr[i][j]->SetTextColor(plot::kMpGrey);
break;
default:
break;
}
t[i][j]->Draw();
}
ctmp->Print(fmt::format("{}{}.png", output_prefix, VarName[i]).c_str());
delete ctmp;
}
// TODO we're not actually getting the resolutions yet
// error for the test result
Q2_resolution_test.error(-1);
//============================== pass the test results ==============================
for(int i = 0 ; i < 4 ; i++){
double width = myf[i][1]->GetParameter(1);
if(myFitPtr[i][1]->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
eic::util::write_test(Q2_resolution_test, fmt::format("{}invar.json", output_prefix));
// That's all!
eic::util::write_test(Tests, fmt::format("{}invar.json", output_prefix));
return 0;
}
Loading