Skip to content
Snippets Groups Projects
Commit fc5a58c0 authored by Ziyue Zhang's avatar Ziyue Zhang
Browse files

WIP: Add t in dvmp; add a place holder function with limited track selection... 

WIP: Add t in dvmp; add a place holder function with limited track selection for invariant quantities in REC
parent f8eba5d4
No related branches found
No related tags found
No related merge requests found
This commit is part of merge request !25. Comments created here will be created in the context of that merge request.
Hide whitespace changes
Inline Side-by-side
benchmarks/dvmp/analysis/dvmp.h
+ 67
− 8
Edit View file @ fc5a58c0
...@@ -32,20 +32,79 @@ namespace util { ...@@ -32,20 +32,79 @@ namespace util {
}; };
// for simu // for simu
inline inv_quant calc_inv_quant_simu(const std::vector<ROOT::Math::PxPyPzMVector>& parts) 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 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 Q2 = - q.Dot(q);
double t = Delta.Dot(Delta);
inv_quant quantities = {nu, Q2, Q2/2./P.mass()/nu, t};
return quantities;
}
// for rec
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;
double nu = q.Dot(P) / P.mass(); // go through all particle combinatorics, calculate the invariant mass
double Q2 = -q.Dot(q); // for each combination, and remember which combination is the closest
inv_quant quantities = {nu, Q2, Q2 / 2. / P.mass() / nu}; // 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; 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; }
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 // for tracking, add later
......
benchmarks/dvmp/analysis/vm_invar.cxx
+ 127
− 68
Edit View file @ fc5a58c0
...@@ -73,6 +73,10 @@ int vm_invar(const std::string& config_name) ...@@ -73,6 +73,10 @@ int vm_invar(const std::string& config_name)
auto momenta_from_tracking = [decay_mass](const std::vector<eic::TrackParametersData>& tracks) { auto momenta_from_tracking = [decay_mass](const std::vector<eic::TrackParametersData>& tracks) {
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);
};
//==================================================================== //====================================================================
...@@ -81,102 +85,157 @@ int vm_invar(const std::string& config_name) ...@@ -81,102 +85,157 @@ 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_simu, {"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
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_rec = d_im.Histo1D({"h_t_rec", ";t;#", 100, -1., 0.}, "t_rec");
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");
auto h_t_sim = d_im.Histo1D({"h_t_sim" , ";t;#", 100, -1., 0.}, "t_sim");
// 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.
{ {
// Print canvas to output file // Print canvas to output file
TCanvas c{"canvas2", "canvas2", 1200, 1200};
TCanvas c{"canvas2", "canvas2", 1800, 600}; c.Divide(2, 2, 0.0001, 0.0001);
c.Divide(3, 1, 0.0001, 0.0001); //pad 1 nu
// pad 1 nu
c.cd(1); c.cd(1);
// gPad->SetLogx(false); //gPad->SetLogx(false);
// gPad->SetLogy(false); //gPad->SetLogy(false);
auto& hnu = *h_nu_sim; auto& hnu_rec = *h_nu_rec;
auto& hnu_sim = *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_rec.GetXaxis()->CenterTitle();
// hnu.GetXaxis()->SetTitle("#times1000"); //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, vm_name, "#nu");
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("reconstructed");
// 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); //gPad->SetLogx(false);
// gPad->SetLogy(false); //gPad->SetLogy(false);
auto& hQ2 = *h_Q2_sim; auto& hQ2_rec = *h_Q2_rec;
auto& hQ2_sim = *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_rec.GetXaxis()->CenterTitle();
//hnu.GetXaxis()->SetTitle("#times1000");
// 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, vm_name, "Q^{2}");
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("reconstructed");
// 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); //gPad->SetLogx(false);
// gPad->SetLogy(false); //gPad->SetLogy(false);
auto& hx = *h_x_sim; auto& hx_rec = *h_x_rec;
auto& hx_sim = *h_x_sim;
// histogram style
hx_rec.SetLineColor(plot::kMpOrange);
hx_rec.SetLineWidth(2);
hx_sim.SetLineColor(plot::kMpBlue);
hx_sim.SetLineWidth(2);
// axes
hx_rec.GetXaxis()->CenterTitle();
//hnu.GetXaxis()->SetTitle("#times1000");
// draw everything
hx_sim.DrawClone("hist");
hx_rec.DrawClone("hist same");
// FIXME hardcoded beam configuration
plot::draw_label(10, 100, detector, vm_name, "x");
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("reconstructed");
tptr3->SetTextColor(plot::kMpOrange);
t3->Draw();
//pad 4 t
c.cd(4);
//gPad->SetLogx(false);
//gPad->SetLogy(false);
auto& ht_rec = *h_t_rec;
auto& ht_sim = *h_t_sim;
// histogram style // histogram style
hx.SetLineColor(plot::kMpBlue); ht_rec.SetLineColor(plot::kMpOrange);
hx.SetLineWidth(2); ht_rec.SetLineWidth(2);
ht_sim.SetLineColor(plot::kMpBlue);
ht_sim.SetLineWidth(2);
// axes // axes
hx.GetXaxis()->CenterTitle(); ht_rec.GetXaxis()->CenterTitle();
//hnu.GetXaxis()->SetTitle("#times1000");
// draw everything // draw everything
hx.DrawClone("hist"); ht_sim.DrawClone("hist");
ht_rec.DrawClone("hist same");
// FIXME hardcoded beam configuration // FIXME hardcoded beam configuration
plot::draw_label(10, 100, detector); plot::draw_label(10, 100, detector, vm_name, "t");
TText* tptr23; TText* tptr4;
auto t23 = new TPaveText(.6, .8417, .9, .925, "NB NDC"); auto t4 = new TPaveText(.6, .8417, .9, .925, "NB NDC");
t23->SetFillColorAlpha(kWhite, 0); t4->SetFillColorAlpha(kWhite, 0);
t23->SetTextFont(43); t4->SetTextFont(43);
t23->SetTextSize(25); t4->SetTextSize(25);
tptr23 = t23->AddText("simulated"); tptr4 = t4->AddText("simulated");
tptr23->SetTextColor(plot::kMpBlue); tptr4->SetTextColor(plot::kMpBlue);
// tptr1 = t1->AddText("reconstructed"); tptr4 = t4->AddText("reconstructed");
// tptr1->SetTextColor(plot::kMpOrange); tptr4->SetTextColor(plot::kMpOrange);
t23->Draw(); t4->Draw();
c.Print(fmt::format("{}InvariantQuantities.png", output_prefix).c_str()); c.Print(fmt::format("{}InvariantQuantities.png", output_prefix).c_str());
} }
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment