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

benchmarks/dvmp/analysis/dvmp.h

@@ -46,7 +46,63 @@ namespace util {
   }
   
   //for tracking
-  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){
+ //   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) {
+ //     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;
+ // }
+ 
+ //for tracking
+  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;
@@ -55,7 +111,9 @@ namespace util {
     // 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.001) continue;
       for (int j = i + 1; j < parts.size(); ++j) {
+        if( fabs(parts[j].mass() - daughter_mass)/daughter_mass > 0.001) continue;
         const double new_mass{(parts[i] + parts[j]).mass()};
         if (fabs(new_mass - pdg_mass) < fabs(best_mass - pdg_mass)) {
           first = i;
@@ -65,7 +123,7 @@ namespace util {
       }
     }
     if (first < 0 || parts.size() < 3 ){
-      inv_quant quantities = {-999., -999., -999., -999.};
+      inv_quant quantities = {-99999., -99999., -99999., -99999.};
       return quantities;
     }
     ROOT::Math::PxPyPzMVector pair_4p(parts[first] + parts[second]);
@@ -100,7 +158,6 @@ namespace util {
     //inv_quant quantities = {-999., -999., -999., -999.};
     return quantities;
   }
-