WIP: Test on new smear

benchmarks/dvmp/analysis/dvmp.h

@@ -79,7 +79,8 @@ namespace util {
     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);
+      //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);
     }
     
@@ -98,10 +99,12 @@ namespace util {
       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);
+        //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);
+        //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))) {
@@ -123,7 +126,8 @@ namespace util {
       //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
-        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);
+        //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;
       //}
@@ -206,7 +210,7 @@ namespace util {
     //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[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();