Skip to content
Snippets Groups Projects
Commit 29052b17 authored by Ziyue Zhang's avatar Ziyue Zhang Committed by Whitney Armstrong
Browse files

Add t in dvmp; add a place holder for InvQuant in REC

parent d35832cf
No related branches found
No related tags found
1 merge request!25Add t in dvmp; add a place holder for InvQuant in REC
Hide whitespace changes
Inline Side-by-side
benchmarks/dvmp/analysis/dvmp.h
+ 70
− 6
View file @ 29052b17
...@@ -28,24 +28,88 @@ namespace util { ...@@ -28,24 +28,88 @@ namespace util {
// for structure functions // for structure functions
struct inv_quant { // add more when needed struct inv_quant { // add more when needed
double nu, Q2, x; double nu, Q2, x, t;
}; };
// for simu // for simu
inline inv_quant calc_inv_quant_simu(const std::vector<ROOT::Math::PxPyPzMVector>& parts) inline inv_quant calc_inv_quant_sim(const std::vector<ROOT::Math::PxPyPzMVector>& parts)
{ {
ROOT::Math::PxPyPzMVector q(parts[0] - parts[2]); ROOT::Math::PxPyPzMVector q(parts[0] - parts[2]);
ROOT::Math::PxPyPzMVector P(parts[3]); ROOT::Math::PxPyPzMVector P(parts[3]);
ROOT::Math::PxPyPzMVector Delta(parts[6] - parts[3]);
double nu = q.Dot(P) / P.mass(); double nu = q.Dot(P) / P.mass();
double Q2 = -q.Dot(q); double Q2 = -q.Dot(q);
inv_quant quantities = {nu, Q2, Q2 / 2. / P.mass() / nu}; double t = Delta.Dot(Delta);
inv_quant quantities = {nu, Q2, Q2 / 2. / P.mass() / nu, t};
return quantities; return quantities;
} }
//for tracking
inline inv_quant calc_inv_quant_rec(const std::vector<ROOT::Math::PxPyPzMVector>& parts, const double pdg_mass){
int first = -1;
int second = -1;
double best_mass = -1;
inline double get_nu_simu(inv_quant quantities) { return quantities.nu / 1000.; } // go through all particle combinatorics, calculate the invariant mass
inline double get_Q2_simu(inv_quant quantities) { return quantities.Q2; } // for each combination, and remember which combination is the closest
inline double get_x_simu(inv_quant quantities) { return quantities.x; } // to the desired pdg_mass
for (int i = 0; i < parts.size(); ++i) {
for (int j = i + 1; j < parts.size(); ++j) {
const double new_mass{(parts[i] + parts[j]).mass()};
if (fabs(new_mass - pdg_mass) < fabs(best_mass - pdg_mass)) {
first = i;
second = j;
best_mass = new_mass;
}
}
}
if (first < 0 || parts.size() < 3 ){
inv_quant quantities = {-999., -999., -999., -999.};
return quantities;
}
ROOT::Math::PxPyPzMVector pair_4p(parts[first] + parts[second]);
ROOT::Math::PxPyPzMVector e1, P;
double e1_Energy = sqrt(10.*10. + get_pdg_mass("electron")*get_pdg_mass("electron"));
double P_Energy = sqrt(100.*100. + get_pdg_mass("proton")*get_pdg_mass("proton"));
e1.SetPxPyPzE(0., 0., -10., e1_Energy);
P.SetPxPyPzE(0., 0., 100., P_Energy);
int scatteredIdx = -1;
float dp = 10.;
for(int i = 0 ; i < parts.size(); i++){
if(i==first || i==second) continue;
ROOT::Math::PxPyPzMVector k_prime(parts[i]);
float ptmp = sqrt(parts[i].px()*parts[i].px() + parts[i].py()*parts[i].py() + parts[i].pz()*parts[i].pz());
if( (k_prime.px()) * (pair_4p.px()) + (k_prime.py()) * (pair_4p.py()) + (k_prime.pz()) * (pair_4p.pz()) > 0. || ptmp >= 10.) continue; //angle between jpsi and scattered electron < pi/2, 3-momentum mag < 10.
if(dp > 10.- ptmp){ //if there are more than one candidate of scattered electron, choose the one with highest 3-momentum mag
scatteredIdx = i;
dp = 10. - ptmp;
}
}
if(scatteredIdx ==-1){
inv_quant quantities = {-999., -999., -999., -999.};
return quantities;
}
ROOT::Math::PxPyPzMVector q(e1 - parts[scatteredIdx]);
ROOT::Math::PxPyPzMVector Delta(q - pair_4p);
double nu = q.Dot(P) / P.mass();
double Q2 = - q.Dot(q);
double t = Delta.Dot(Delta);
inv_quant quantities = {nu, Q2, Q2/2./P.mass()/nu, t};
//inv_quant quantities = {-999., -999., -999., -999.};
return quantities;
}
inline double get_nu(inv_quant quantities) { return quantities.nu / 1000.; }
inline double get_Q2(inv_quant quantities) { return quantities.Q2; }
inline double get_x(inv_quant quantities) { return quantities.x; }
inline double get_t(inv_quant quantities) { return quantities.t; }
// for tracking, add later // for tracking, add later
......
benchmarks/dvmp/analysis/vm_invar.cxx
+ 111
− 61
View file @ 29052b17
...@@ -74,6 +74,10 @@ int vm_invar(const std::string& config_name) ...@@ -74,6 +74,10 @@ int vm_invar(const std::string& config_name)
return util::momenta_from_tracking(tracks, decay_mass); return util::momenta_from_tracking(tracks, decay_mass);
}; };
auto calc_inv_quant_rec = [vm_mass](const std::vector<ROOT::Math::PxPyPzMVector>& parts) {
return util::calc_inv_quant_rec(parts, vm_mass);
};
//==================================================================== //====================================================================
// Define analysis flow // Define analysis flow
...@@ -81,10 +85,16 @@ int vm_invar(const std::string& config_name) ...@@ -81,10 +85,16 @@ int vm_invar(const std::string& config_name)
.Define("N", "p_rec.size()") .Define("N", "p_rec.size()")
.Define("p_sim", util::momenta_from_simulation, {"mcparticles2"}) .Define("p_sim", util::momenta_from_simulation, {"mcparticles2"})
//================================================================ //================================================================
.Define("invariant_quantities", util::calc_inv_quant_simu, {"p_sim"}) .Define("invariant_quantities_rec", calc_inv_quant_rec, {"p_rec"})
.Define("nu_sim", util::get_nu_simu, {"invariant_quantities"}) .Define("invariant_quantities_sim", util::calc_inv_quant_sim, {"p_sim"})
.Define("Q2_sim", util::get_Q2_simu, {"invariant_quantities"}) .Define("nu_rec", util::get_nu, {"invariant_quantities_rec"})
.Define("x_sim", util::get_x_simu, {"invariant_quantities"}); .Define("Q2_rec", util::get_Q2, {"invariant_quantities_rec"})
.Define("x_rec", util::get_x, {"invariant_quantities_rec"})
.Define("t_rec", util::get_t, {"invariant_quantities_rec"})
.Define("nu_sim", util::get_nu, {"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 // Define output histograms
...@@ -92,6 +102,12 @@ int vm_invar(const std::string& config_name) ...@@ -92,6 +102,12 @@ int vm_invar(const std::string& config_name)
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};#", 100, 0., 15.}, "Q2_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"); auto h_x_sim = d_im.Histo1D({"h_x_sim", ";x;#", 100, 0., 0.1}, "x_sim");
auto h_t_sim = d_im.Histo1D({"h_t_sim", ";t;#", 100, -1., 0.}, "t_sim");
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};#", 100, 0., 15.}, "Q2_rec");
auto h_x_rec = d_im.Histo1D({"h_x_rec", ";x;#", 100, 0., 0.1}, "x_rec");
auto h_t_rec = d_im.Histo1D({"h_t_rec", ";t;#", 100, -1., 0.}, "t_rec");
// 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
...@@ -100,85 +116,119 @@ int vm_invar(const std::string& config_name) ...@@ -100,85 +116,119 @@ int vm_invar(const std::string& config_name)
// Print canvas to output file // Print canvas to output file
TCanvas c{"canvas2", "canvas2", 1800, 600}; TCanvas c{"canvas2", "canvas2", 1200, 900};
c.Divide(3, 1, 0.0001, 0.0001); c.Divide(2, 2, 0.0001, 0.0001);
// pad 1 nu // pad 1 nu
c.cd(1); c.cd(1);
// gPad->SetLogx(false); auto& hnu_rec = *h_nu_rec;
// gPad->SetLogy(false); auto& hnu_sim = *h_nu_sim;
auto& hnu = *h_nu_sim;
// histogram style // histogram style
hnu.SetLineColor(plot::kMpBlue); hnu_rec.SetLineColor(plot::kMpOrange);
hnu.SetLineWidth(2); hnu_rec.SetLineWidth(2);
hnu_sim.SetLineColor(plot::kMpBlue);
hnu_sim.SetLineWidth(2);
// axes // axes
hnu.GetXaxis()->CenterTitle(); hnu_sim.GetXaxis()->CenterTitle();
// hnu.GetXaxis()->SetTitle("#times1000");
// draw everything // draw everything
hnu.DrawClone("hist"); hnu_sim.DrawClone("hist");
hnu_rec.DrawClone("hist same");
// FIXME hardcoded beam configuration // FIXME hardcoded beam configuration
plot::draw_label(10, 100, detector); plot::draw_label(10, 100, detector);
TText* tptr21; TText* tptr1;
auto t21 = new TPaveText(.6, .8417, .9, .925, "NB NDC"); auto t1 = new TPaveText(.6, .8417, .9, .925, "NB NDC");
t21->SetFillColorAlpha(kWhite, 0); t1->SetFillColorAlpha(kWhite, 0);
t21->SetTextFont(43); t1->SetTextFont(43);
t21->SetTextSize(25); t1->SetTextSize(25);
tptr21 = t21->AddText("simulated"); tptr1 = t1->AddText("simulated");
tptr21->SetTextColor(plot::kMpBlue); tptr1->SetTextColor(plot::kMpBlue);
// tptr1 = t1->AddText("reconstructed"); tptr1 = t1->AddText("rec(PlaceHolder)");
// tptr1->SetTextColor(plot::kMpOrange); tptr1->SetTextColor(plot::kMpOrange);
t21->Draw(); t1->Draw();
// pad 2 Q2 // pad 2 Q2
c.cd(2); c.cd(2);
// gPad->SetLogx(false); auto& hQ2_rec = *h_Q2_rec;
// gPad->SetLogy(false); auto& hQ2_sim = *h_Q2_sim;
auto& hQ2 = *h_Q2_sim;
// histogram style // histogram style
hQ2.SetLineColor(plot::kMpBlue); hQ2_rec.SetLineColor(plot::kMpOrange);
hQ2.SetLineWidth(2); hQ2_rec.SetLineWidth(2);
hQ2_sim.SetLineColor(plot::kMpBlue);
hQ2_sim.SetLineWidth(2);
// axes // axes
hQ2.GetXaxis()->CenterTitle(); hQ2_sim.GetXaxis()->CenterTitle();
// draw everything // draw everything
hQ2.DrawClone("hist"); hQ2_sim.DrawClone("hist");
hQ2_rec.DrawClone("hist same");
// FIXME hardcoded beam configuration // FIXME hardcoded beam configuration
plot::draw_label(10, 100, detector); plot::draw_label(10, 100, detector);
TText* tptr22; TText* tptr2;
auto t22 = new TPaveText(.6, .8417, .9, .925, "NB NDC"); auto t2 = new TPaveText(.6, .8417, .9, .925, "NB NDC");
t22->SetFillColorAlpha(kWhite, 0); t2->SetFillColorAlpha(kWhite, 0);
t22->SetTextFont(43); t2->SetTextFont(43);
t22->SetTextSize(25); t2->SetTextSize(25);
tptr22 = t22->AddText("simulated"); tptr2 = t2->AddText("simulated");
tptr22->SetTextColor(plot::kMpBlue); tptr2->SetTextColor(plot::kMpBlue);
// tptr1 = t1->AddText("reconstructed"); tptr2 = t2->AddText("rec(PlaceHolder)");
// tptr1->SetTextColor(plot::kMpOrange); tptr2->SetTextColor(plot::kMpOrange);
t22->Draw(); t2->Draw();
// pad 1 nu // pad 3 x
c.cd(3); c.cd(3);
// gPad->SetLogx(false); auto& hx_rec = *h_x_rec;
// gPad->SetLogy(false); auto& hx_sim = *h_x_sim;
auto& hx = *h_x_sim;
// histogram style // histogram style
hx.SetLineColor(plot::kMpBlue); hx_rec.SetLineColor(plot::kMpOrange);
hx.SetLineWidth(2); hx_rec.SetLineWidth(2);
hx_sim.SetLineColor(plot::kMpBlue);
hx_sim.SetLineWidth(2);
// axes // axes
hx.GetXaxis()->CenterTitle(); hx_sim.GetXaxis()->CenterTitle();
// draw everything // draw everything
hx.DrawClone("hist"); hx_sim.DrawClone("hist");
hx_rec.DrawClone("hist same");
// FIXME hardcoded beam configuration // FIXME hardcoded beam configuration
plot::draw_label(10, 100, detector); plot::draw_label(10, 100, detector);
TText* tptr23; TText* tptr3;
auto t23 = new TPaveText(.6, .8417, .9, .925, "NB NDC"); auto t3 = new TPaveText(.6, .8417, .9, .925, "NB NDC");
t23->SetFillColorAlpha(kWhite, 0); t3->SetFillColorAlpha(kWhite, 0);
t23->SetTextFont(43); t3->SetTextFont(43);
t23->SetTextSize(25); t3->SetTextSize(25);
tptr23 = t23->AddText("simulated"); tptr3 = t3->AddText("simulated");
tptr23->SetTextColor(plot::kMpBlue); tptr3->SetTextColor(plot::kMpBlue);
// tptr1 = t1->AddText("reconstructed"); tptr3 = t3->AddText("rec(PlaceHolder)");
// tptr1->SetTextColor(plot::kMpOrange); tptr3->SetTextColor(plot::kMpOrange);
t23->Draw(); t3->Draw();
// pad 3 x
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(2);
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();
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 // TODO we're not actually getting the resolutions yet
......
include/util.h
+ 2
− 0
View file @ 29052b17
...@@ -48,6 +48,8 @@ namespace util { ...@@ -48,6 +48,8 @@ namespace util {
return 3.0969; return 3.0969;
} else if (part == "upsilon") { } else if (part == "upsilon") {
return 9.49630; return 9.49630;
} else if (part == "proton"){
return 0.938272;
} else { } else {
throw unknown_particle_error{part}; throw unknown_particle_error{part};
} }
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment