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

Add J/psi resolution; re-write RC with Dummy instead of rc tracks

parent 17ebacf1
No related branches found
No related tags found
1 merge request!27Add J/psi resolution; re-write RC with Dummy instead of rc tracks
Hide whitespace changes
Inline Side-by-side
benchmarks/dvmp/analysis/dvmp.h
+ 32
− 33
View file @ 27efeb40
...@@ -17,13 +17,6 @@ ...@@ -17,13 +17,6 @@
// promoted to the top-level util library // promoted to the top-level util library
namespace util { namespace util {
// Calculate the 4-vector sum of a given pair of particles
inline ROOT::Math::PxPyPzMVector
get_sum(const std::pair<ROOT::Math::PxPyPzMVector, ROOT::Math::PxPyPzMVector>& particle_pair)
{
return (particle_pair.first + particle_pair.second);
}
//======================================================================================================== //========================================================================================================
// for structure functions // for structure functions
...@@ -36,26 +29,29 @@ namespace util { ...@@ -36,26 +29,29 @@ namespace util {
{ {
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]); //ROOT::Math::PxPyPzMVector Delta(parts[6] - parts[3]);//exact jpsi
ROOT::Math::PxPyPzMVector Delta(parts[0] - parts[2] - parts[7] - parts[8]);//jpsi->l l' gamma, ignore gamma
double nu = q.Dot(P) / P.mass(); double nu = q.Dot(P) / P.mass();
double Q2 = -q.Dot(q); double Q2 = -q.Dot(q);
double t = Delta.Dot(Delta); double t = Delta.Dot(Delta);
inv_quant quantities = {nu, Q2, Q2 / 2. / P.mass() / nu, t}; inv_quant quantities = {nu, Q2, Q2/2./P.mass()/nu, t};
return quantities; return quantities;
} }
//for tracking //for Dummy rc
inline inv_quant calc_inv_quant_rec(const std::vector<ROOT::Math::PxPyPzMVector>& parts, const double pdg_mass){ inline inv_quant calc_inv_quant_rec(const std::vector<ROOT::Math::PxPyPzMVector>& parts, const double pdg_mass, const double daughter_mass){
int first = -1; int first = -1;
int second = -1; int second = -1;
double best_mass = -1; double best_mass = -1;
// go through all particle combinatorics, calculate the invariant mass // go through all particle combinatorics, calculate the invariant mass
// for each combination, and remember which combination is the closest // for each combination, and remember which combination is the closest
// to the desired pdg_mass // to the desired pdg_mass
for (int i = 0; i < parts.size(); ++i) { for (int i = 0; i < parts.size(); ++i) {
if( fabs(parts[i].mass() - daughter_mass)/daughter_mass > 0.01) continue;
for (int j = i + 1; j < parts.size(); ++j) { for (int j = i + 1; j < parts.size(); ++j) {
if( fabs(parts[j].mass() - daughter_mass)/daughter_mass > 0.01) continue;
const double new_mass{(parts[i] + parts[j]).mass()}; const double new_mass{(parts[i] + parts[j]).mass()};
if (fabs(new_mass - pdg_mass) < fabs(best_mass - pdg_mass)) { if (fabs(new_mass - pdg_mass) < fabs(best_mass - pdg_mass)) {
first = i; first = i;
...@@ -64,30 +60,39 @@ namespace util { ...@@ -64,30 +60,39 @@ namespace util {
} }
} }
} }
if (first < 0 || parts.size() < 3 ){
inv_quant quantities = {-999., -999., -999., -999.};
if (first < 0 || parts.size() < 3 ){ //fewer than 3 candidates
inv_quant quantities = {-99999., -99999., -99999., -99999.};
return quantities; return quantities;
} }
//construct the beam kinematics
ROOT::Math::PxPyPzMVector pair_4p(parts[first] + parts[second]); ROOT::Math::PxPyPzMVector pair_4p(parts[first] + parts[second]);
ROOT::Math::PxPyPzMVector e1, P; ROOT::Math::PxPyPzMVector e1, P;
double e1_Energy = sqrt(10.*10. + get_pdg_mass("electron")*get_pdg_mass("electron")); //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")); //double P_Energy = sqrt(100.*100. + get_pdg_mass("proton")*get_pdg_mass("proton"));
e1.SetPxPyPzE(0., 0., -10., e1_Energy); double e1_Energy = sqrt((1.305e-8 - 10.)*(1.305e-8 - 10.) + (0.0005109+9.888e-8)*(0.0005109+9.888e-8));
P.SetPxPyPzE(0., 0., 100., P_Energy); double P_Energy = sqrt((99.995598 + 1.313e-7)*(99.995598 + 1.313e-7) + (0.938272-1.23e-12)*(0.938272-1.23e-12));
e1.SetPxPyPzE(0., 0., 1.305e-8 - 10., e1_Energy);
P.SetPxPyPzE(0., 0., 99.995598 + 1.313e-7, P_Energy);
int scatteredIdx = -1; int scatteredIdx = -1;
float dp = 10.;
//float dp = 10.;
for(int i = 0 ; i < parts.size(); i++){ for(int i = 0 ; i < parts.size(); i++){
if(i==first || i==second) continue; if(i==first || i==second) continue; //skip the paired leptons
ROOT::Math::PxPyPzMVector k_prime(parts[i]); //if( fabs(parts[i].mass() - get_pdg_mass("electron"))/get_pdg_mass("electron") > 0.01) continue;
float ptmp = sqrt(parts[i].px()*parts[i].px() + parts[i].py()*parts[i].py() + parts[i].pz()*parts[i].pz()); if( fabs(parts[i].mass() - 0.0005109 - 9.888e-8)/(0.0005109 + 9.888e-8) > 0.01) continue;
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. //ROOT::Math::PxPyPzMVector k_prime(parts[i]); //scattered
if(dp > 10.- ptmp){ //if there are more than one candidate of scattered electron, choose the one with highest 3-momentum mag //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; scatteredIdx = i;
dp = 10. - ptmp; // dp = 10. - ptmp;
} //}
} }
if(scatteredIdx ==-1){ if(scatteredIdx ==-1){
inv_quant quantities = {-999., -999., -999., -999.}; inv_quant quantities = {-99999., -99999., -99999., -99999.};
return quantities; return quantities;
} }
ROOT::Math::PxPyPzMVector q(e1 - parts[scatteredIdx]); ROOT::Math::PxPyPzMVector q(e1 - parts[scatteredIdx]);
...@@ -97,15 +102,9 @@ namespace util { ...@@ -97,15 +102,9 @@ namespace util {
double Q2 = - q.Dot(q); double Q2 = - q.Dot(q);
double t = Delta.Dot(Delta); double t = Delta.Dot(Delta);
inv_quant quantities = {nu, Q2, Q2/2./P.mass()/nu, t}; 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_nu(inv_quant quantities) { return quantities.nu / 1000.; }
inline double get_Q2(inv_quant quantities) { return quantities.Q2; } inline double get_Q2(inv_quant quantities) { return quantities.Q2; }
inline double get_x(inv_quant quantities) { return quantities.x; } inline double get_x(inv_quant quantities) { return quantities.x; }
......
benchmarks/dvmp/analysis/vm_invar.cxx
+ 8
− 11
View file @ 27efeb40
...@@ -70,18 +70,15 @@ int vm_invar(const std::string& config_name) ...@@ -70,18 +70,15 @@ int vm_invar(const std::string& config_name)
// utility lambda functions to bind the vector meson and decay particle // utility lambda functions to bind the vector meson and decay particle
// types // types
auto momenta_from_tracking = [decay_mass](const std::vector<eic::TrackParametersData>& tracks) {
return util::momenta_from_tracking(tracks, decay_mass); auto calc_inv_quant_rec = [vm_mass, decay_mass](const std::vector<ROOT::Math::PxPyPzMVector>& parts) {
}; return util::calc_inv_quant_rec(parts, vm_mass, 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
auto d_im = d.Define("p_rec", momenta_from_tracking, {"outputTrackParameters"}) auto d_im = d.Define("p_rec", util::momenta_RC, {"DummyReconstructedParticles"})
.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"})
//================================================================ //================================================================
...@@ -124,7 +121,7 @@ int vm_invar(const std::string& config_name) ...@@ -124,7 +121,7 @@ int vm_invar(const std::string& config_name)
auto& hnu_sim = *h_nu_sim; auto& hnu_sim = *h_nu_sim;
// histogram style // histogram style
hnu_rec.SetLineColor(plot::kMpOrange); hnu_rec.SetLineColor(plot::kMpOrange);
hnu_rec.SetLineWidth(2); hnu_rec.SetLineWidth(1);
hnu_sim.SetLineColor(plot::kMpBlue); hnu_sim.SetLineColor(plot::kMpBlue);
hnu_sim.SetLineWidth(2); hnu_sim.SetLineWidth(2);
// axes // axes
...@@ -151,7 +148,7 @@ int vm_invar(const std::string& config_name) ...@@ -151,7 +148,7 @@ int vm_invar(const std::string& config_name)
auto& hQ2_sim = *h_Q2_sim; auto& hQ2_sim = *h_Q2_sim;
// histogram style // histogram style
hQ2_rec.SetLineColor(plot::kMpOrange); hQ2_rec.SetLineColor(plot::kMpOrange);
hQ2_rec.SetLineWidth(2); hQ2_rec.SetLineWidth(1);
hQ2_sim.SetLineColor(plot::kMpBlue); hQ2_sim.SetLineColor(plot::kMpBlue);
hQ2_sim.SetLineWidth(2); hQ2_sim.SetLineWidth(2);
// axes // axes
...@@ -178,7 +175,7 @@ int vm_invar(const std::string& config_name) ...@@ -178,7 +175,7 @@ int vm_invar(const std::string& config_name)
auto& hx_sim = *h_x_sim; auto& hx_sim = *h_x_sim;
// histogram style // histogram style
hx_rec.SetLineColor(plot::kMpOrange); hx_rec.SetLineColor(plot::kMpOrange);
hx_rec.SetLineWidth(2); hx_rec.SetLineWidth(1);
hx_sim.SetLineColor(plot::kMpBlue); hx_sim.SetLineColor(plot::kMpBlue);
hx_sim.SetLineWidth(2); hx_sim.SetLineWidth(2);
// axes // axes
...@@ -205,7 +202,7 @@ int vm_invar(const std::string& config_name) ...@@ -205,7 +202,7 @@ int vm_invar(const std::string& config_name)
auto& ht_sim = *h_t_sim; auto& ht_sim = *h_t_sim;
// histogram style // histogram style
ht_rec.SetLineColor(plot::kMpOrange); ht_rec.SetLineColor(plot::kMpOrange);
ht_rec.SetLineWidth(2); ht_rec.SetLineWidth(1);
ht_sim.SetLineColor(plot::kMpBlue); ht_sim.SetLineColor(plot::kMpBlue);
ht_sim.SetLineWidth(2); ht_sim.SetLineWidth(2);
// axes // axes
......
benchmarks/dvmp/analysis/vm_mass.cxx
+ 43
− 21
View file @ 27efeb40
...@@ -14,6 +14,8 @@ ...@@ -14,6 +14,8 @@
#include <nlohmann/json.hpp> #include <nlohmann/json.hpp>
#include <string> #include <string>
#include <vector> #include <vector>
#include "eicd/ReconstructedParticleCollection.h"
#include "eicd/ReconstructedParticleData.h"
// Run VM invariant-mass-based benchmarks on an input reconstruction file for // 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) // a desired vector meson (e.g. jpsi) and a desired decay particle (e.g. muon)
...@@ -22,6 +24,19 @@ ...@@ -22,6 +24,19 @@
// TODO: I think it would be better to pass small json configuration file to // 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. // the test, instead of this ever-expanding list of function arguments.
// FIXME: MC does not trace back into particle history. Need to fix that // FIXME: MC does not trace back into particle history. Need to fix that
//double RBW(double*x, double*par){
// double mean = par[0];
// double GAMMA = par[1];
// double N = par[2];
// double gamma = mean*TMath::Sqrt(mean*mean + GAMMA*GAMMA);
// double k = 2.*mean*GAMMA*gamma/TMath::Pi()*TMath::Sqrt(2./(mean*mean + gamma));
// double eval = N*k/((x[0]*x[0]-mean*mean)*(x[0]*x[0]-mean*mean) + mean*mean*GAMMA*GAMMA);
// return(eval);
//}
//double fFlat(double*x, double*par){
// return(par[0]);
//}
int vm_mass(const std::string& config_name) int vm_mass(const std::string& config_name)
{ {
// read our configuration // read our configuration
...@@ -72,23 +87,20 @@ int vm_mass(const std::string& config_name) ...@@ -72,23 +87,20 @@ int vm_mass(const std::string& config_name)
// utility lambda functions to bind the vector meson and decay particle // utility lambda functions to bind the vector meson and decay particle
// types // types
auto momenta_from_tracking = [decay_mass](const std::vector<eic::TrackParametersData>& tracks) {
return util::momenta_from_tracking(tracks, decay_mass); auto find_decay_pair = [vm_mass, decay_mass](const std::vector<ROOT::Math::PxPyPzMVector>& parts) {
}; return util::find_decay_pair(parts, vm_mass, decay_mass);
auto find_decay_pair = [vm_mass](const std::vector<ROOT::Math::PxPyPzMVector>& parts) {
return util::find_decay_pair(parts, vm_mass);
}; };
// util::PrintGeant4(mcparticles2); // util::PrintGeant4(mcparticles2);
// Define analysis flow // Define analysis flow
auto d_im = d.Define("p_rec", momenta_from_tracking, {"outputTrackParameters"}) auto d_im = d.Define("p_rec", util::momenta_RC, {"DummyReconstructedParticles"}) //using dummy rc
.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("decay_pair_rec", find_decay_pair, {"p_rec"}) .Define("decay_pair_rec", find_decay_pair, {"p_rec"})
.Define("decay_pair_sim", find_decay_pair, {"p_sim"}) .Define("decay_pair_sim", find_decay_pair, {"p_sim"})
.Define("p_vm_rec", "decay_pair_rec.first + decay_pair_rec.second") .Define("p_vm_rec", "decay_pair_rec.first + decay_pair_rec.second")
.Define("p_vm_sim", "decay_pair_sim.first + decay_pair_sim.second") .Define("p_vm_sim", "decay_pair_sim.first + decay_pair_sim.second")
//.Define("p_vm_sim", util::get_sum, {"decay_pair_sim"})
.Define("mass_rec", "p_vm_rec.M()") .Define("mass_rec", "p_vm_rec.M()")
.Define("mass_sim", "p_vm_sim.M()") .Define("mass_sim", "p_vm_sim.M()")
.Define("pt_rec", "p_vm_rec.pt()") .Define("pt_rec", "p_vm_rec.pt()")
...@@ -99,19 +111,18 @@ int vm_mass(const std::string& config_name) ...@@ -99,19 +111,18 @@ int vm_mass(const std::string& config_name)
.Define("eta_sim", "p_vm_sim.eta()"); .Define("eta_sim", "p_vm_sim.eta()");
// Define output histograms // Define output histograms
auto h_im_rec = //auto h_im_rec = d_im.Histo1D({"h_im_rec", ";m_{ll'} (GeV/c^{2});#", (int)(vm_mass+0.5)*2*100, 0., 2.*(int)(vm_mass+0.5)}, "mass_rec"); //real rec
d_im.Histo1D({"h_im_rec", ";m_{ll'} (GeV/c^{2});#", 100, -1.1, vm_mass + 5}, "mass_rec"); auto h_im_rec = d_im.Histo1D({"h_im_rec", ";m_{ll'} (GeV/c^{2});#", 50, 1., 5.}, "mass_rec");//for dummy_rec
auto h_im_sim = auto h_im_sim = d_im.Histo1D({"h_im_sim", ";m_{ll'} (GeV/c^{2});#", 50, 1., 5.}, "mass_sim");
d_im.Histo1D({"h_im_sim", ";m_{ll'} (GeV/c^{2});#", 100, -1.1, vm_mass + 5}, "mass_sim");
auto h_pt_rec = d_im.Histo1D({"h_pt_rec", ";p_{T} (GeV/c);#", 400, 0., 40.}, "pt_rec"); auto h_pt_rec = d_im.Histo1D({"h_pt_rec", ";p_{T} (GeV/c);#", 50, 0., 10.}, "pt_rec");
auto h_pt_sim = d_im.Histo1D({"h_pt_sim", ";p_{T} (GeV/c);#", 400, 0., 40.}, "pt_sim"); auto h_pt_sim = d_im.Histo1D({"h_pt_sim", ";p_{T} (GeV/c);#", 50, 0., 10.}, "pt_sim");
auto h_phi_rec = d_im.Histo1D({"h_phi_rec", ";#phi_{ll'};#", 90, -M_PI, M_PI}, "phi_rec"); auto h_phi_rec = d_im.Histo1D({"h_phi_rec", ";#phi_{ll'};#", 45, -M_PI, M_PI}, "phi_rec");
auto h_phi_sim = d_im.Histo1D({"h_phi_sim", ";#phi_{ll'};#", 90, -M_PI, M_PI}, "phi_sim"); auto h_phi_sim = d_im.Histo1D({"h_phi_sim", ";#phi_{ll'};#", 45, -M_PI, M_PI}, "phi_sim");
auto h_eta_rec = d_im.Histo1D({"h_eta_rec", ";#eta_{ll'};#", 1000, -5., 5.}, "eta_rec"); auto h_eta_rec = d_im.Histo1D({"h_eta_rec", ";#eta_{ll'};#", 50, -2., 2.}, "eta_rec");
auto h_eta_sim = d_im.Histo1D({"h_eta_sim", ";#eta_{ll'};#", 1000, -5., 5.}, "eta_sim"); auto h_eta_sim = d_im.Histo1D({"h_eta_sim", ";#eta_{ll'};#", 50, -2., 2.}, "eta_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
...@@ -129,7 +140,7 @@ int vm_mass(const std::string& config_name) ...@@ -129,7 +140,7 @@ int vm_mass(const std::string& config_name)
h11.SetLineColor(plot::kMpBlue); h11.SetLineColor(plot::kMpBlue);
h11.SetLineWidth(2); h11.SetLineWidth(2);
h12.SetLineColor(plot::kMpOrange); h12.SetLineColor(plot::kMpOrange);
h12.SetLineWidth(2); h12.SetLineWidth(1);
// axes // axes
h11.GetXaxis()->CenterTitle(); h11.GetXaxis()->CenterTitle();
h11.GetYaxis()->CenterTitle(); h11.GetYaxis()->CenterTitle();
...@@ -159,13 +170,24 @@ int vm_mass(const std::string& config_name) ...@@ -159,13 +170,24 @@ int vm_mass(const std::string& config_name)
h21.SetLineColor(plot::kMpBlue); h21.SetLineColor(plot::kMpBlue);
h21.SetLineWidth(2); h21.SetLineWidth(2);
h22.SetLineColor(plot::kMpOrange); h22.SetLineColor(plot::kMpOrange);
h22.SetLineWidth(2); h22.SetLineWidth(1);
// axes // axes
h21.GetXaxis()->CenterTitle(); h21.GetXaxis()->CenterTitle();
h21.GetYaxis()->CenterTitle(); h21.GetYaxis()->CenterTitle();
// draw everything // draw everything
TF1* mfPt = new TF1("mfPt", "[0]*exp(-[1]*x)", 1.2, 10.);
mfPt->SetParameters(5., 1.);
mfPt->SetParLimits(0, 0., 1000000.);
mfPt->SetParLimits(1, 0., 1000000.);
mfPt->SetNpx(1000);
mfPt->SetLineColor(2);
mfPt->SetLineStyle(7);
h21.Fit(mfPt, "", "", 1.2, 10.);
h21.DrawClone("hist"); h21.DrawClone("hist");
h22.DrawClone("hist same"); h22.DrawClone("hist same");
mfPt->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;
...@@ -189,7 +211,7 @@ int vm_mass(const std::string& config_name) ...@@ -189,7 +211,7 @@ int vm_mass(const std::string& config_name)
h31.SetLineColor(plot::kMpBlue); h31.SetLineColor(plot::kMpBlue);
h31.SetLineWidth(2); h31.SetLineWidth(2);
h32.SetLineColor(plot::kMpOrange); h32.SetLineColor(plot::kMpOrange);
h32.SetLineWidth(2); h32.SetLineWidth(1);
// axes // axes
h31.GetXaxis()->CenterTitle(); h31.GetXaxis()->CenterTitle();
h31.GetYaxis()->CenterTitle(); h31.GetYaxis()->CenterTitle();
...@@ -219,7 +241,7 @@ int vm_mass(const std::string& config_name) ...@@ -219,7 +241,7 @@ int vm_mass(const std::string& config_name)
h41.SetLineColor(plot::kMpBlue); h41.SetLineColor(plot::kMpBlue);
h41.SetLineWidth(2); h41.SetLineWidth(2);
h42.SetLineColor(plot::kMpOrange); h42.SetLineColor(plot::kMpOrange);
h42.SetLineWidth(2); h42.SetLineWidth(1);
// axes // axes
h41.GetXaxis()->CenterTitle(); h41.GetXaxis()->CenterTitle();
h41.GetYaxis()->CenterTitle(); h41.GetYaxis()->CenterTitle();
......
include/util.h
+ 18
− 3
View file @ 27efeb40
...@@ -10,11 +10,14 @@ ...@@ -10,11 +10,14 @@
#include <limits> #include <limits>
#include <string> #include <string>
#include <vector> #include <vector>
#include "TF1.h"
#include <Math/Vector4D.h> #include <Math/Vector4D.h>
#include "dd4pod/Geant4ParticleCollection.h" #include "dd4pod/Geant4ParticleCollection.h"
#include "eicd/TrackParametersCollection.h" #include "eicd/TrackParametersCollection.h"
#include "eicd/ReconstructedParticleCollection.h"
#include "eicd/ReconstructedParticleData.h"
namespace util { namespace util {
...@@ -56,7 +59,7 @@ namespace util { ...@@ -56,7 +59,7 @@ namespace util {
} }
// Get a vector of 4-momenta from raw tracking info, using an externally // Get a vector of 4-momenta from raw tracking info, using an externally
// provided particle mass assumption. // provided particle mass assumption. //outputTrackParameters
inline auto momenta_from_tracking(const std::vector<eic::TrackParametersData>& tracks, inline auto momenta_from_tracking(const std::vector<eic::TrackParametersData>& tracks,
const double mass) const double mass)
{ {
...@@ -75,7 +78,17 @@ namespace util { ...@@ -75,7 +78,17 @@ namespace util {
}); });
return momenta; return momenta;
} }
//=====================================================================================
inline auto momenta_RC(const std::vector<eic::ReconstructedParticleData>& parts)
{
std::vector<ROOT::Math::PxPyPzMVector> momenta{parts.size()};
// transform our raw tracker info into proper 4-momenta
std::transform(parts.begin(), parts.end(), momenta.begin(), [](const auto& part) {
return ROOT::Math::PxPyPzMVector{part.p.x, part.p.y, part.p.z, part.mass};
});
return momenta;
}
//=====================================================================================
// Get a vector of 4-momenta from the simulation data. // Get a vector of 4-momenta from the simulation data.
// TODO: Add PID selector (maybe using ranges?) // TODO: Add PID selector (maybe using ranges?)
inline auto momenta_from_simulation(const std::vector<dd4pod::Geant4ParticleData>& parts) inline auto momenta_from_simulation(const std::vector<dd4pod::Geant4ParticleData>& parts)
...@@ -91,7 +104,7 @@ namespace util { ...@@ -91,7 +104,7 @@ namespace util {
// Find the decay pair candidates from a vector of particles (parts), // Find the decay pair candidates from a vector of particles (parts),
// with invariant mass closest to a desired value (pdg_mass) // with invariant mass closest to a desired value (pdg_mass)
inline std::pair<ROOT::Math::PxPyPzMVector, ROOT::Math::PxPyPzMVector> inline std::pair<ROOT::Math::PxPyPzMVector, ROOT::Math::PxPyPzMVector>
find_decay_pair(const std::vector<ROOT::Math::PxPyPzMVector>& parts, const double pdg_mass) find_decay_pair(const std::vector<ROOT::Math::PxPyPzMVector>& parts, const double pdg_mass, const double daughter_mass)
{ {
int first = -1; int first = -1;
int second = -1; int second = -1;
...@@ -101,7 +114,9 @@ namespace util { ...@@ -101,7 +114,9 @@ namespace util {
// for each combination, and remember which combination is the closest // for each combination, and remember which combination is the closest
// to the desired pdg_mass // to the desired pdg_mass
for (size_t i = 0; i < parts.size(); ++i) { for (size_t i = 0; i < parts.size(); ++i) {
if( fabs(parts[i].mass() - daughter_mass)/daughter_mass > 0.01) continue;
for (size_t j = i + 1; j < parts.size(); ++j) { for (size_t j = i + 1; j < parts.size(); ++j) {
if( fabs(parts[j].mass() - daughter_mass)/daughter_mass > 0.01) continue;
const double new_mass{(parts[i] + parts[j]).mass()}; const double new_mass{(parts[i] + parts[j]).mass()};
if (fabs(new_mass - pdg_mass) < fabs(best_mass - pdg_mass)) { if (fabs(new_mass - pdg_mass) < fabs(best_mass - pdg_mass)) {
first = i; first = i;
......
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