diff --git a/benchmarks/dvmp/analysis/dvmp.h b/benchmarks/dvmp/analysis/dvmp.h
index 3baf03ad4d9d8213d52ada928b34ed74631fa0e4..20e6b9537a5a7c4478993f1bc74c1a3fb92b042e 100644
--- a/benchmarks/dvmp/analysis/dvmp.h
+++ b/benchmarks/dvmp/analysis/dvmp.h
@@ -17,13 +17,6 @@
// promoted to the top-level util library
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
@@ -36,26 +29,29 @@ namespace util {
{
ROOT::Math::PxPyPzMVector q(parts[0] - parts[2]);
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 Q2 = -q.Dot(q);
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;
}
- //for tracking
- inline inv_quant calc_inv_quant_rec(const std::vector<ROOT::Math::PxPyPzMVector>& parts, const double pdg_mass){
+ //for Dummy rc
+ 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 second = -1;
double best_mass = -1;
-
+
// go through all particle combinatorics, calculate the invariant mass
// for each combination, and remember which combination is the closest
// to the desired pdg_mass
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) {
+ if( fabs(parts[j].mass() - daughter_mass)/daughter_mass > 0.01) continue;
const double new_mass{(parts[i] + parts[j]).mass()};
if (fabs(new_mass - pdg_mass) < fabs(best_mass - pdg_mass)) {
first = i;
@@ -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;
}
+
+ //construct the beam kinematics
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);
+ //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 e1_Energy = sqrt((1.305e-8 - 10.)*(1.305e-8 - 10.) + (0.0005109+9.888e-8)*(0.0005109+9.888e-8));
+ 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;
- float dp = 10.;
+
+ //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
+ if(i==first || i==second) continue; //skip the paired leptons
+ //if( fabs(parts[i].mass() - get_pdg_mass("electron"))/get_pdg_mass("electron") > 0.01) continue;
+ if( fabs(parts[i].mass() - 0.0005109 - 9.888e-8)/(0.0005109 + 9.888e-8) > 0.01) continue;
+ //ROOT::Math::PxPyPzMVector k_prime(parts[i]); //scattered
+ //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;
- }
+ // dp = 10. - ptmp;
+ //}
}
if(scatteredIdx ==-1){
- inv_quant quantities = {-999., -999., -999., -999.};
+ inv_quant quantities = {-99999., -99999., -99999., -99999.};
return quantities;
}
ROOT::Math::PxPyPzMVector q(e1 - parts[scatteredIdx]);
@@ -97,15 +102,9 @@ namespace util {
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; }
diff --git a/benchmarks/dvmp/analysis/vm_invar.cxx b/benchmarks/dvmp/analysis/vm_invar.cxx
index edaa03ad2adb1546f15a1d547c5e75c9685d3e46..21be3710106ac4495b03c26ae5a9d2eed4499b80 100644
--- a/benchmarks/dvmp/analysis/vm_invar.cxx
+++ b/benchmarks/dvmp/analysis/vm_invar.cxx
@@ -70,18 +70,15 @@ int vm_invar(const std::string& config_name)
// utility lambda functions to bind the vector meson and decay particle
// 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](const std::vector<ROOT::Math::PxPyPzMVector>& parts) {
- return util::calc_inv_quant_rec(parts, vm_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);
};
//====================================================================
// 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("p_sim", util::momenta_from_simulation, {"mcparticles2"})
//================================================================
@@ -124,7 +121,7 @@ int vm_invar(const std::string& config_name)
auto& hnu_sim = *h_nu_sim;
// histogram style
hnu_rec.SetLineColor(plot::kMpOrange);
- hnu_rec.SetLineWidth(2);
+ hnu_rec.SetLineWidth(1);
hnu_sim.SetLineColor(plot::kMpBlue);
hnu_sim.SetLineWidth(2);
// axes
@@ -151,7 +148,7 @@ int vm_invar(const std::string& config_name)
auto& hQ2_sim = *h_Q2_sim;
// histogram style
hQ2_rec.SetLineColor(plot::kMpOrange);
- hQ2_rec.SetLineWidth(2);
+ hQ2_rec.SetLineWidth(1);
hQ2_sim.SetLineColor(plot::kMpBlue);
hQ2_sim.SetLineWidth(2);
// axes
@@ -178,7 +175,7 @@ int vm_invar(const std::string& config_name)
auto& hx_sim = *h_x_sim;
// histogram style
hx_rec.SetLineColor(plot::kMpOrange);
- hx_rec.SetLineWidth(2);
+ hx_rec.SetLineWidth(1);
hx_sim.SetLineColor(plot::kMpBlue);
hx_sim.SetLineWidth(2);
// axes
@@ -205,7 +202,7 @@ int vm_invar(const std::string& config_name)
auto& ht_sim = *h_t_sim;
// histogram style
ht_rec.SetLineColor(plot::kMpOrange);
- ht_rec.SetLineWidth(2);
+ ht_rec.SetLineWidth(1);
ht_sim.SetLineColor(plot::kMpBlue);
ht_sim.SetLineWidth(2);
// axes
diff --git a/benchmarks/dvmp/analysis/vm_mass.cxx b/benchmarks/dvmp/analysis/vm_mass.cxx
index 37c5d386386049b517de0ccc1a5d457936f3ac9c..48b05db406d0853205e66d72333a7e2257c53a0e 100644
--- a/benchmarks/dvmp/analysis/vm_mass.cxx
+++ b/benchmarks/dvmp/analysis/vm_mass.cxx
@@ -14,6 +14,8 @@
#include <nlohmann/json.hpp>
#include <string>
#include <vector>
+#include "eicd/ReconstructedParticleCollection.h"
+#include "eicd/ReconstructedParticleData.h"
// 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)
@@ -22,6 +24,19 @@
// 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.
// 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)
{
// read our configuration
@@ -72,23 +87,20 @@ int vm_mass(const std::string& config_name)
// utility lambda functions to bind the vector meson and decay particle
// 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](const std::vector<ROOT::Math::PxPyPzMVector>& parts) {
- return util::find_decay_pair(parts, vm_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);
};
// util::PrintGeant4(mcparticles2);
// 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("p_sim", util::momenta_from_simulation, {"mcparticles2"})
.Define("decay_pair_rec", find_decay_pair, {"p_rec"})
.Define("decay_pair_sim", find_decay_pair, {"p_sim"})
.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", util::get_sum, {"decay_pair_sim"})
.Define("mass_rec", "p_vm_rec.M()")
.Define("mass_sim", "p_vm_sim.M()")
.Define("pt_rec", "p_vm_rec.pt()")
@@ -99,19 +111,18 @@ int vm_mass(const std::string& config_name)
.Define("eta_sim", "p_vm_sim.eta()");
// Define output histograms
- auto h_im_rec =
- d_im.Histo1D({"h_im_rec", ";m_{ll'} (GeV/c^{2});#", 100, -1.1, vm_mass + 5}, "mass_rec");
- auto h_im_sim =
- d_im.Histo1D({"h_im_sim", ";m_{ll'} (GeV/c^{2});#", 100, -1.1, vm_mass + 5}, "mass_sim");
+ //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
+ 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 = d_im.Histo1D({"h_im_sim", ";m_{ll'} (GeV/c^{2});#", 50, 1., 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_sim = d_im.Histo1D({"h_pt_sim", ";p_{T} (GeV/c);#", 400, 0., 40.}, "pt_sim");
+ 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);#", 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_sim = d_im.Histo1D({"h_phi_sim", ";#phi_{ll'};#", 90, -M_PI, M_PI}, "phi_sim");
+ 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'};#", 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_sim = d_im.Histo1D({"h_eta_sim", ";#eta_{ll'};#", 1000, -5., 5.}, "eta_sim");
+ 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'};#", 50, -2., 2.}, "eta_sim");
// Plot our histograms.
// TODO: to start I'm explicitly plotting the histograms, but want to
@@ -129,7 +140,7 @@ int vm_mass(const std::string& config_name)
h11.SetLineColor(plot::kMpBlue);
h11.SetLineWidth(2);
h12.SetLineColor(plot::kMpOrange);
- h12.SetLineWidth(2);
+ h12.SetLineWidth(1);
// axes
h11.GetXaxis()->CenterTitle();
h11.GetYaxis()->CenterTitle();
@@ -159,13 +170,24 @@ int vm_mass(const std::string& config_name)
h21.SetLineColor(plot::kMpBlue);
h21.SetLineWidth(2);
h22.SetLineColor(plot::kMpOrange);
- h22.SetLineWidth(2);
+ h22.SetLineWidth(1);
// axes
h21.GetXaxis()->CenterTitle();
h21.GetYaxis()->CenterTitle();
// 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");
h22.DrawClone("hist same");
+ mfPt->Draw("same");
+
+
// FIXME hardcoded beam configuration
plot::draw_label(10, 100, detector);
TText* tptr2;
@@ -189,7 +211,7 @@ int vm_mass(const std::string& config_name)
h31.SetLineColor(plot::kMpBlue);
h31.SetLineWidth(2);
h32.SetLineColor(plot::kMpOrange);
- h32.SetLineWidth(2);
+ h32.SetLineWidth(1);
// axes
h31.GetXaxis()->CenterTitle();
h31.GetYaxis()->CenterTitle();
@@ -219,7 +241,7 @@ int vm_mass(const std::string& config_name)
h41.SetLineColor(plot::kMpBlue);
h41.SetLineWidth(2);
h42.SetLineColor(plot::kMpOrange);
- h42.SetLineWidth(2);
+ h42.SetLineWidth(1);
// axes
h41.GetXaxis()->CenterTitle();
h41.GetYaxis()->CenterTitle();
diff --git a/include/util.h b/include/util.h
index 6a24b293c26963b4999fc63df821ab3418c694a8..eddf13933a08e4d0c3b0de952b421fb1904e8ea8 100644
--- a/include/util.h
+++ b/include/util.h
@@ -10,11 +10,14 @@
#include <limits>
#include <string>
#include <vector>
+#include "TF1.h"
#include <Math/Vector4D.h>
#include "dd4pod/Geant4ParticleCollection.h"
#include "eicd/TrackParametersCollection.h"
+#include "eicd/ReconstructedParticleCollection.h"
+#include "eicd/ReconstructedParticleData.h"
namespace util {
@@ -56,7 +59,7 @@ namespace util {
}
// 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,
const double mass)
{
@@ -75,7 +78,17 @@ namespace util {
});
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.
// TODO: Add PID selector (maybe using ranges?)
inline auto momenta_from_simulation(const std::vector<dd4pod::Geant4ParticleData>& parts)
@@ -91,7 +104,7 @@ namespace util {
// Find the decay pair candidates from a vector of particles (parts),
// with invariant mass closest to a desired value (pdg_mass)
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 second = -1;
@@ -101,7 +114,9 @@ namespace util {
// for each combination, and remember which combination is the closest
// to the desired pdg_mass
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) {
+ if( fabs(parts[j].mass() - daughter_mass)/daughter_mass > 0.01) continue;
const double new_mass{(parts[i] + parts[j]).mass()};
if (fabs(new_mass - pdg_mass) < fabs(best_mass - pdg_mass)) {
first = i;