diff --git a/benchmarks/dvmp/analysis/dvmp.h b/benchmarks/dvmp/analysis/dvmp.h
index 371b454f8d92f4c0cd5b17edd4a8f5e1ecc77e77..350467d8e8229b591b8757eeec14317bd4d9d27d 100644
--- a/benchmarks/dvmp/analysis/dvmp.h
+++ b/benchmarks/dvmp/analysis/dvmp.h
@@ -32,20 +32,79 @@ namespace util {
   };
 
   // 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 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();
-    double    Q2         = -q.Dot(q);
-    inv_quant quantities = {nu, Q2, Q2 / 2. / P.mass() / nu};
+    // 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) {
+      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;          }
 
-  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
 
diff --git a/benchmarks/dvmp/analysis/vm_invar.cxx b/benchmarks/dvmp/analysis/vm_invar.cxx
index bbc4db2739073fb604c4bbeb7b2d5d0746d03b3a..14dd695286e86f124f9c476c1b920eb5726858f1 100644
--- a/benchmarks/dvmp/analysis/vm_invar.cxx
+++ b/benchmarks/dvmp/analysis/vm_invar.cxx
@@ -73,6 +73,10 @@ int vm_invar(const std::string& config_name)
   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);
+  };
 
   //====================================================================
 
@@ -81,102 +85,157 @@ int vm_invar(const std::string& config_name)
                   .Define("N", "p_rec.size()")
                   .Define("p_sim", util::momenta_from_simulation, {"mcparticles2"})
                   //================================================================
-                  .Define("invariant_quantities", util::calc_inv_quant_simu, {"p_sim"})
-                  .Define("nu_sim", util::get_nu_simu, {"invariant_quantities"})
-                  .Define("Q2_sim", util::get_Q2_simu, {"invariant_quantities"})
-                  .Define("x_sim", util::get_x_simu, {"invariant_quantities"});
-  //================================================================
+                  .Define("invariant_quantities_rec", calc_inv_quant_rec, {"p_rec"})
+                  .Define("invariant_quantities_sim", util::calc_inv_quant_simu, {"p_sim"})
+                  .Define("nu_rec" , util::get_nu, {"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("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_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.
   // TODO: to start I'm explicitly plotting the histograms, but want to
   // factorize out the plotting code moving forward.
   {
 
     // Print canvas to output file
-
-    TCanvas c{"canvas2", "canvas2", 1800, 600};
-    c.Divide(3, 1, 0.0001, 0.0001);
-    // pad 1 nu
+TCanvas c{"canvas2", "canvas2", 1200, 1200};
+    c.Divide(2, 2, 0.0001, 0.0001);
+    //pad 1 nu
     c.cd(1);
-    // gPad->SetLogx(false);
-    // gPad->SetLogy(false);
-    auto& hnu = *h_nu_sim;
+    //gPad->SetLogx(false);
+    //gPad->SetLogy(false);
+    auto& hnu_rec = *h_nu_rec;
+    auto& hnu_sim = *h_nu_sim;
     // histogram style
-    hnu.SetLineColor(plot::kMpBlue);
-    hnu.SetLineWidth(2);
+    hnu_rec.SetLineColor(plot::kMpOrange);
+    hnu_rec.SetLineWidth(2);
+    hnu_sim.SetLineColor(plot::kMpBlue);
+    hnu_sim.SetLineWidth(2);
     // axes
-    hnu.GetXaxis()->CenterTitle();
-    // hnu.GetXaxis()->SetTitle("#times1000");
+    hnu_rec.GetXaxis()->CenterTitle();
+    //hnu.GetXaxis()->SetTitle("#times1000");
     // draw everything
-    hnu.DrawClone("hist");
+    hnu_sim.DrawClone("hist");
+    hnu_rec.DrawClone("hist same");
     // FIXME hardcoded beam configuration
-    plot::draw_label(10, 100, detector);
-    TText* tptr21;
-    auto   t21 = new TPaveText(.6, .8417, .9, .925, "NB NDC");
-    t21->SetFillColorAlpha(kWhite, 0);
-    t21->SetTextFont(43);
-    t21->SetTextSize(25);
-    tptr21 = t21->AddText("simulated");
-    tptr21->SetTextColor(plot::kMpBlue);
-    // tptr1 = t1->AddText("reconstructed");
-    // tptr1->SetTextColor(plot::kMpOrange);
-    t21->Draw();
-
-    // pad 2 Q2
+    plot::draw_label(10, 100, detector, vm_name, "#nu");
+    TText* tptr1;
+    auto t1 = new TPaveText(.6, .8417, .9, .925, "NB NDC");
+    t1->SetFillColorAlpha(kWhite, 0);
+    t1->SetTextFont(43);
+    t1->SetTextSize(25);
+    tptr1 = t1->AddText("simulated");
+    tptr1->SetTextColor(plot::kMpBlue);
+    tptr1 = t1->AddText("reconstructed");
+    tptr1->SetTextColor(plot::kMpOrange);
+    t1->Draw();
+    
+    //pad 2 Q2
     c.cd(2);
-    // gPad->SetLogx(false);
-    // gPad->SetLogy(false);
-    auto& hQ2 = *h_Q2_sim;
+    //gPad->SetLogx(false);
+    //gPad->SetLogy(false);
+    auto& hQ2_rec = *h_Q2_rec;
+    auto& hQ2_sim = *h_Q2_sim;
     // histogram style
-    hQ2.SetLineColor(plot::kMpBlue);
-    hQ2.SetLineWidth(2);
+    hQ2_rec.SetLineColor(plot::kMpOrange);
+    hQ2_rec.SetLineWidth(2);
+    hQ2_sim.SetLineColor(plot::kMpBlue);
+    hQ2_sim.SetLineWidth(2);
     // axes
-    hQ2.GetXaxis()->CenterTitle();
+    hQ2_rec.GetXaxis()->CenterTitle();
+    //hnu.GetXaxis()->SetTitle("#times1000");
     // draw everything
-    hQ2.DrawClone("hist");
+    hQ2_sim.DrawClone("hist");
+    hQ2_rec.DrawClone("hist same");
     // FIXME hardcoded beam configuration
-    plot::draw_label(10, 100, detector);
-    TText* tptr22;
-    auto   t22 = new TPaveText(.6, .8417, .9, .925, "NB NDC");
-    t22->SetFillColorAlpha(kWhite, 0);
-    t22->SetTextFont(43);
-    t22->SetTextSize(25);
-    tptr22 = t22->AddText("simulated");
-    tptr22->SetTextColor(plot::kMpBlue);
-    // tptr1 = t1->AddText("reconstructed");
-    // tptr1->SetTextColor(plot::kMpOrange);
-    t22->Draw();
-
-    // pad 1 nu
+    plot::draw_label(10, 100, detector, vm_name, "Q^{2}");
+    TText* tptr2;
+    auto t2 = new TPaveText(.6, .8417, .9, .925, "NB NDC");
+    t2->SetFillColorAlpha(kWhite, 0);
+    t2->SetTextFont(43);
+    t2->SetTextSize(25);
+    tptr2 = t2->AddText("simulated");
+    tptr2->SetTextColor(plot::kMpBlue);
+    tptr2 = t2->AddText("reconstructed");
+    tptr2->SetTextColor(plot::kMpOrange);
+    t2->Draw();
+
+    
+    //pad 3 x
     c.cd(3);
-    // gPad->SetLogx(false);
-    // gPad->SetLogy(false);
-    auto& hx = *h_x_sim;
+    //gPad->SetLogx(false);
+    //gPad->SetLogy(false);
+    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
-    hx.SetLineColor(plot::kMpBlue);
-    hx.SetLineWidth(2);
+    ht_rec.SetLineColor(plot::kMpOrange);
+    ht_rec.SetLineWidth(2);
+    ht_sim.SetLineColor(plot::kMpBlue);
+    ht_sim.SetLineWidth(2);
     // axes
-    hx.GetXaxis()->CenterTitle();
+    ht_rec.GetXaxis()->CenterTitle();
+    //hnu.GetXaxis()->SetTitle("#times1000");
     // draw everything
-    hx.DrawClone("hist");
+    ht_sim.DrawClone("hist");
+    ht_rec.DrawClone("hist same");
     // FIXME hardcoded beam configuration
-    plot::draw_label(10, 100, detector);
-    TText* tptr23;
-    auto   t23 = new TPaveText(.6, .8417, .9, .925, "NB NDC");
-    t23->SetFillColorAlpha(kWhite, 0);
-    t23->SetTextFont(43);
-    t23->SetTextSize(25);
-    tptr23 = t23->AddText("simulated");
-    tptr23->SetTextColor(plot::kMpBlue);
-    // tptr1 = t1->AddText("reconstructed");
-    // tptr1->SetTextColor(plot::kMpOrange);
-    t23->Draw();
+    plot::draw_label(10, 100, detector, vm_name, "t");
+    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("reconstructed");
+    tptr4->SetTextColor(plot::kMpOrange);
+    t4->Draw();
 
     c.Print(fmt::format("{}InvariantQuantities.png", output_prefix).c_str());
   }