diff --git a/benchmarks/dvmp/analysis/dvmp.h b/benchmarks/dvmp/analysis/dvmp.h
index 20e6b9537a5a7c4478993f1bc74c1a3fb92b042e..40f3e82b7689cd010eaf930ca8464c9bfc642390 100644
--- a/benchmarks/dvmp/analysis/dvmp.h
+++ b/benchmarks/dvmp/analysis/dvmp.h
@@ -19,95 +19,134 @@ namespace util {
//========================================================================================================
// for structure functions
-
struct inv_quant { // add more when needed
- double nu, Q2, x, t;
+ double nu, Q2, x, y, t;
};
-
- // for simu
- 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 P(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};
- return quantities;
+ //import Geant4 and set the wanted particles in the intended order
+ //0:e0 1:p0 2:e1 3:p1 4:recoil system (without p1) 5:l1 from 4 6:l2 from 4
+ inline auto momenta_sort_sim(const std::vector<dd4pod::Geant4ParticleData>& parts, std::string_view mother, std::string_view daughter){//mother and daughter are not used yet; will be useful when generater is different and/or when the mcparticles doesn't follow the same order in all events
+ std::vector<ROOT::Math::PxPyPzMVector> momenta{7};
+ int order_map[7] = {0, 3, 2, 6, 5, 7, 8};
+ for(int i = 0 ; i < 7 ; i++){
+ double px = parts[order_map[i]].psx;
+ double py = parts[order_map[i]].psy;
+ double pz = parts[order_map[i]].psz;
+ double mass = parts[order_map[i]].mass;
+ double e = sqrt(px*px + py*py + pz*pz + mass*mass);
+ momenta[i].SetPxPyPzE(px, py, pz, e);
+ }
+ //for(int i = 0 ; i < 7 ; i++){
+ // cout<<Form("sim, idx = %d, P(px, py, pz, mass) = (%f, %f, %f, %f)", i, momenta[i].px(), momenta[i].py(), momenta[i].pz(), momenta[i].mass())<<endl;
+ //}
+ //cout<<"========================================="<<endl;
+ return momenta;
}
- //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){
+ //import Reconstructed particles and set the wanted particles in the intended order========================
+ inline auto momenta_sort_rec(const std::vector<eic::ReconstructedParticleData>& parts, std::string_view mother, std::string_view daughter){
+ std::vector<ROOT::Math::PxPyPzMVector> momenta{7};
+ //0:e0 1:p0 2:e1 3:p1 4:recoil system (without p1) 5:l1 from recoil decay 6:l2 from recoil decay
+ for(int i = 0 ; i < 7 ; i++) momenta[i].SetPxPyPzE(0., 0., 0., 0.); //initialize as all 0
+
+ //manually set incoming electron and proton;
+ double e0_mass = get_pdg_mass("electron");
+ double e0_pz = 1.305e-8 - 10.;
+ momenta[0].SetPxPyPzE(0., 0., e0_pz, sqrt(e0_mass*e0_mass + e0_pz*e0_pz));
+ double p0_mass = get_pdg_mass("proton");
+ double p0_pz = 99.995598 + 1.313e-7 + 8.783e-11;
+ momenta[1].SetPxPyPzE(0., 0., p0_pz, sqrt(p0_mass*p0_mass + p0_pz*p0_pz));
+
+ //FIXME search for recoil proton, add feature to deal with multiple protons in one event
+ for(int i = 0 ; i < parts.size() ; i++){
+ if(parts[i].pid!=2212) continue;
+ //px, py, pz, e, mass are not consistent for smeared dummy rc. this is a temp solution, replacing non-smeared energy by the re-calculated energy
+ //double energy_tmp = sqrt(parts[i].p.x*parts[i].p.x + parts[i].p.y*parts[i].p.y + parts[i].p.z*parts[i].p.z + parts[i].mass*parts[i].mass);//tmpsolution
+ double energy_tmp = parts[i].energy;
+ momenta[3].SetPxPyPzE(parts[i].p.x, parts[i].p.y, parts[i].p.z, energy_tmp);
+ }
+
+ //search for di-lepton pair for the decay in recoil
+ int daughter_pid = -1; //unsigned
+ if(daughter == "electron"){
+ daughter_pid = 11;
+ }else if(daughter == "muon"){
+ daughter_pid = 13;
+ }
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)) {
+ double best_mass = -999.;
+ for (int i = 0 ; i < parts.size() ; ++i) {
+ if( fabs(parts[i].pid)!=daughter_pid) continue;
+ for (int j = i + 1 ; j < parts.size() ; ++j) {
+ if( parts[j].pid!= - parts[i].pid) continue;
+ ROOT::Math::PxPyPzMVector lpt_1(0.,0.,0.,0.);
+ //double energy_tmp1 = sqrt(parts[i].p.x*parts[i].p.x + parts[i].p.y*parts[i].p.y + parts[i].p.z*parts[i].p.z + parts[i].mass*parts[i].mass);//tmpsolution
+ double energy_tmp1 = parts[i].energy;
+ lpt_1.SetPxPyPzE(parts[i].p.x, parts[i].p.y, parts[i].p.z, energy_tmp1);
+ ROOT::Math::PxPyPzMVector lpt_2(0.,0.,0.,0.);
+ //double energy_tmp2 = sqrt(parts[j].p.x*parts[j].p.x + parts[j].p.y*parts[j].p.y + parts[j].p.z*parts[j].p.z + parts[j].mass*parts[j].mass);//tmpsolution
+ double energy_tmp2 = parts[j].energy;
+ lpt_2.SetPxPyPzE(parts[j].p.x, parts[j].p.y, parts[j].p.z, energy_tmp2);
+ const double new_mass{(lpt_1 + lpt_2).mass()};
+ if (fabs(new_mass - get_pdg_mass(mother)) < fabs(best_mass - get_pdg_mass(mother))) {
first = i;
second = j;
best_mass = new_mass;
+ momenta[5].SetPxPyPzE(parts[i].p.x, parts[i].p.y, parts[i].p.z, energy_tmp1);
+ momenta[6].SetPxPyPzE(parts[j].p.x, parts[j].p.y, parts[j].p.z, energy_tmp2);
}
+
}
}
-
- 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"));
- 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;
-
+ //FIXME search for scattered electron, need improvement with more complex events
//float dp = 10.;
for(int i = 0 ; i < parts.size(); i++){
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(parts[i].pid != 11) continue;
+ //float ptmp = sqrt(parts[i].p.x*parts[i].p.x + parts[i].p.y*parts[i].p.y + parts[i].p.z*parts[i].p.z);
//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;
+ //double energy_tmp = sqrt(parts[i].p.x*parts[i].p.x + parts[i].p.y*parts[i].p.y + parts[i].p.z*parts[i].p.z + parts[i].mass*parts[i].mass);//tmpsolution
+ double energy_tmp = parts[i].energy;
+ momenta[2].SetPxPyPzE(parts[i].p.x, parts[i].p.y, parts[i].p.z, energy_tmp);
// dp = 10. - ptmp;
//}
}
- if(scatteredIdx ==-1){
- inv_quant quantities = {-99999., -99999., -99999., -99999.};
- 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};
+ //for(int i = 0 ; i < 7 ; i++){
+ // cout<<Form("dum, idx = %d, P(px, py, pz, mass) = (%f, %f, %f, %f)", i, momenta[i].px(), momenta[i].py(), momenta[i].pz(), momenta[i].mass())<<endl;
+ //}
+ //cout<<"========================================="<<endl;
+ return momenta;
+ }
+
+ inline inv_quant calc_inv_quant(const std::vector<ROOT::Math::PxPyPzMVector>& parts)
+ {
+ //0:e0 1:p0 2:e1 3:p1 4:recoil system (without p1) 5:l1 from 4 6:l2 from 4
+ ROOT::Math::PxPyPzMVector q(parts[0] - parts[2]);
+ ROOT::Math::PxPyPzMVector k(parts[0]);
+ ROOT::Math::PxPyPzMVector P(parts[1]);
+ ROOT::Math::PxPyPzMVector Delta(parts[3] - parts[1]);//exact
+ ROOT::Math::PxPyPzMVector Delta_prime(parts[0] - parts[2] - parts[5] - parts[6]);//exclude gamma radiation in jpsi decay
+
+ double nu = q.Dot(P) / P.mass();
+ double Q2 = -q.Dot(q);
+ double t = 0.;
+ //if(parts[4].px() == 0. && parts[4].py() == 0. && parts[4].pz() == 0. && parts[4].mass() == 0.){
+ //t = Delta_prime.Dot(Delta_prime);
+ //}else{
+ t = Delta.Dot(Delta);
+ //}
+ double y = q.Dot(P) / k.Dot(P);
+ inv_quant quantities = {nu, Q2, Q2/2./P.mass()/nu, y, t};
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_y(inv_quant quantities) { return quantities.y; }
inline double get_t(inv_quant quantities) { return quantities.t; }
// for tracking, add later
diff --git a/benchmarks/dvmp/analysis/vm_invar.cxx b/benchmarks/dvmp/analysis/vm_invar.cxx
index 21be3710106ac4495b03c26ae5a9d2eed4499b80..5c544e2be49d801889135bd67be94eb86120439e 100644
--- a/benchmarks/dvmp/analysis/vm_invar.cxx
+++ b/benchmarks/dvmp/analysis/vm_invar.cxx
@@ -71,41 +71,46 @@ int vm_invar(const std::string& config_name)
// utility lambda functions to bind the vector meson and decay particle
// types
- 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 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
- auto d_im = d.Define("p_rec", util::momenta_RC, {"DummyReconstructedParticles"})
- .Define("N", "p_rec.size()")
- .Define("p_sim", util::momenta_from_simulation, {"mcparticles2"})
- //================================================================
- .Define("invariant_quantities_rec", calc_inv_quant_rec, {"p_rec"})
- .Define("invariant_quantities_sim", util::calc_inv_quant_sim, {"p_sim"})
- .Define("nu_rec", util::get_nu, {"invariant_quantities_rec"})
+ 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()")
+ .Define("invariant_quantities_rec", util::calc_inv_quant, {"p_rec_sorted"})
+ .Define("invariant_quantities_sim", util::calc_inv_quant, {"p_sim_sorted"})
+ .Define("y_rec", util::get_y, {"invariant_quantities_rec"})
.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("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};#", 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");
+ //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");
+
+
+ //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");
- 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.
// TODO: to start I'm explicitly plotting the histograms, but want to
// factorize out the plotting code moving forward.
@@ -117,18 +122,18 @@ int vm_invar(const std::string& config_name)
c.Divide(2, 2, 0.0001, 0.0001);
// pad 1 nu
c.cd(1);
- auto& hnu_rec = *h_nu_rec;
- auto& hnu_sim = *h_nu_sim;
+ auto& hy_rec = *h_y_rec;
+ auto& hy_sim = *h_y_sim;
// histogram style
- hnu_rec.SetLineColor(plot::kMpOrange);
- hnu_rec.SetLineWidth(1);
- hnu_sim.SetLineColor(plot::kMpBlue);
- hnu_sim.SetLineWidth(2);
+ hy_rec.SetLineColor(plot::kMpOrange);
+ hy_rec.SetLineWidth(1);
+ hy_sim.SetLineColor(plot::kMpBlue);
+ hy_sim.SetLineWidth(2);
// axes
- hnu_sim.GetXaxis()->CenterTitle();
+ hy_sim.GetXaxis()->CenterTitle();
// draw everything
- hnu_sim.DrawClone("hist");
- hnu_rec.DrawClone("hist same");
+ hy_sim.DrawClone("hist");
+ hy_rec.DrawClone("hist same");
// FIXME hardcoded beam configuration
plot::draw_label(10, 100, detector);
TText* tptr1;
@@ -196,7 +201,7 @@ int vm_invar(const std::string& config_name)
tptr3->SetTextColor(plot::kMpOrange);
t3->Draw();
- // pad 3 x
+ // pad 4 t
c.cd(4);
auto& ht_rec = *h_t_rec;
auto& ht_sim = *h_t_sim;
diff --git a/benchmarks/dvmp/analysis/vm_mass.cxx b/benchmarks/dvmp/analysis/vm_mass.cxx
index 48b05db406d0853205e66d72333a7e2257c53a0e..9d9b32b1aa8d0b51113483587c9b8d552f373aec 100644
--- a/benchmarks/dvmp/analysis/vm_mass.cxx
+++ b/benchmarks/dvmp/analysis/vm_mass.cxx
@@ -175,17 +175,17 @@ int vm_mass(const std::string& config_name)
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.);
+ //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");
+ //mfPt->Draw("same");
// FIXME hardcoded beam configuration
diff --git a/options/tracker_reconstruction.py b/options/tracker_reconstruction.py
index 50fff6b829cb92e79c0c2d4a68a4c12a68d5755e..6bef5d7b817f3bc0910830e4bb17b7598f3c38bc 100644
--- a/options/tracker_reconstruction.py
+++ b/options/tracker_reconstruction.py
@@ -57,8 +57,7 @@ podioinput = PodioInput("PodioReader",
dummy = MC2DummyParticle("MC2Dummy",
inputCollection="mcparticles",
outputCollection="DummyReconstructedParticles",
- smearing=0.0
- )
+ smearing = 0.0)
## copiers to get around input --> output copy bug. Note the "2" appended to the output collection.
copier = MCCopier("MCCopier",