WIP: rec_ordered(rec) and sim_ordered(sim) debug

benchmarks/dvmp/analysis/dvmp.h

@@ -54,11 +54,6 @@ namespace util {
       double e = sqrt(px*px + py*py + pz*pz + mass*mass);
       momenta[i].SetPxPyPzE(px, py, pz, e);
     }
-    //cout<<Form("sim_ordered, P0 = (%f, %f, %f, %f), P1 = (%f, %f, %f, %f)", momenta[1].px(), momenta[1].py(), momenta[1].pz(), momenta[1].mass(), momenta[3].px(), momenta[3].py(), momenta[3].pz(), momenta[3].mass())<<endl;
-    for(int i = 0 ; i < 7 ; i++){
-        cout<<Form("Sim, %d, (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;
   }
   
@@ -82,11 +77,8 @@ namespace util {
       //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);
       momenta[3].SetPxPyPzE(parts[i].p.x,  parts[i].p.y,  parts[i].p.z,  energy_tmp);
-      //cout<<Form("rec_ordered, p1(px, py, pz, mass) = (%f, %f, %f, %f)", momenta[3].px(), momenta[3].py(), momenta[3].pz(), momenta[3].mass())<<endl;
     }
     
-    //cout<<Form("rec_ordered, P0 = (%f, %f, %f, %f), P1 = (%f, %f, %f, %f)", momenta[1].px(), momenta[1].py(), momenta[1].pz(), momenta[1].mass(), momenta[3].px(), momenta[3].py(), momenta[3].pz(), momenta[3].mass())<<endl;
-    
     //search for di-lepton pair for the decay in recoil
     int daughter_pid = -1;    //unsigned
     if(daughter == "electron"){
@@ -118,10 +110,6 @@ namespace util {
         
       }
     }
-    /*if(first !=  -1){
-      momenta[5].SetPxPyPzE(parts[first].p.x,  parts[first].p.y,  parts[first].p.z,  parts[first].energy);
-      momenta[6].SetPxPyPzE(parts[second].p.x, parts[second].p.y, parts[second].p.z, parts[second].energy);
-    }*/
     
     //FIXME search for scattered electron, need improvement with more complex events
     //float dp = 10.;
@@ -136,12 +124,6 @@ namespace util {
       //  dp = 10. - ptmp;
       //}
     }
-    for(int i = 0 ; i < 7 ; i++){
-        cout<<Form("Dum, %d, (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;
   }
   
@@ -217,11 +199,17 @@ namespace util {
     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
+    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 = Delta.Dot(Delta);
+    double t = 0.;
+    if(parts[4].px() == 0. && parts[4].py() == 0. && parts[4].pz() == 0. && parts[4].mass() == 0.){
+        t = Delta.Dot(Delta);
+    }else{
+        t = Delta_prime.Dot(Delta_prime);
+    }
+    
     inv_quant quantities = {nu, Q2, Q2/2./P.mass()/nu, t};
     return quantities;
   }