WIP: Test on PID instead of pmass

benchmarks/dvmp/analysis/dvmp.h

@@ -120,16 +120,6 @@ namespace util {
     return momenta;
   }
   
-  
-  
-  inline double test_jpsi_pt(const std::vector<ROOT::Math::PxPyPzMVector>& momenta){
-    double px = momenta[4].px();
-    double py = momenta[4].py();
-    double pt = sqrt(px*px + py*py);
-    return pt;
-  }
-  
-  
  //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;
@@ -196,6 +186,22 @@ namespace util {
     return quantities;
   }
   
+  inline inv_quant calc_inv_quant_both(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 P(parts[1]);
+    ROOT::Math::PxPyPzMVector Delta(parts[3] - parts[1]);//jpsi->l l' gamma, ignore gamma
+
+    //ROOT::Math::PxPyPzMVector Delta(parts[6] - parts[3]);//exact jpsi
+    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;
+  }
+  
+  
   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; }

benchmarks/dvmp/analysis/vm_invar.cxx

@@ -84,7 +84,8 @@ int vm_invar(const std::string& config_name)
 
   // Define analysis flow
   auto d_im = d.Define("p_rec", util::momenta_RC, {"DummyReconstructedParticles"})
-                  .Define("test_sim_ordered", momenta_ordered_sim, {"mcparticles2"})
+                  .Define("p_sim_ordered", momenta_ordered_sim, {"mcparticles2"})
+                  .Define("p_rec_ordered", momenta_ordered_rec, {"DummyReconstructedParticles"})
                   .Define("N", "p_rec.size()")
                   .Define("p_sim", util::momenta_from_simulation, {"mcparticles2"})
                   //================================================================
@@ -98,7 +99,8 @@ int vm_invar(const std::string& config_name)
                   .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("test_pt", util::test_jpsi_pt, {"test_sim_ordered"});
+                  .Define("invariant_quantities_rec", calc_inv_quant_both, {"p_rec_ordered"})
+                  .Define("invariant_quantities_sim", calc_inv_quant_both, {"p_sim_ordered"});
   //================================================================
 
   // Define output histograms
@@ -113,8 +115,6 @@ int vm_invar(const std::string& config_name)
   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");
   
-  auto h_test_pt = d_im.Histo1D({"h_test_pt", "h_test_pt", 50, 0., 10.}, "test_pt");
-  
   // Plot our histograms.
   // TODO: to start I'm explicitly plotting the histograms, but want to
   // factorize out the plotting code moving forward.
@@ -205,9 +205,9 @@ 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_rec = *h_t_rec;
     auto& ht_sim = *h_t_sim;
     // histogram style
     ht_rec.SetLineColor(plot::kMpOrange);
@@ -230,25 +230,7 @@ int vm_invar(const std::string& config_name)
     tptr4->SetTextColor(plot::kMpBlue);
     tptr4 = t4->AddText("rec(PlaceHolder)");
     tptr4->SetTextColor(plot::kMpOrange);
-    t4->Draw();*/
-    auto& htest = *h_test_pt;
-    htest.SetLineColor(plot::kMpBlue);
-    htest.SetLineWidth(2);
-    htest.GetXaxis()->CenterTitle();
-    htest.GetYaxis()->CenterTitle();
-    htest.DrawClone("hist");
-    plot::draw_label(10, 100, detector);
-    TText* tptrtest;
-    auto   ttest = new TPaveText(.6, .8417, .9, .925, "NB NDC");
-    ttest->SetFillColorAlpha(kWhite, 0);
-    ttest->SetTextFont(43);
-    ttest->SetTextSize(25);
-    tptrtest = ttest->AddText("simulated");
-    tptrtest->SetTextColor(plot::kMpBlue);
-    ttest->Draw();
-    
-    
-    
+    t4->Draw();
     
     c.Print(fmt::format("{}InvariantQuantities.png", output_prefix).c_str());