From 4ff1c0870c0c01e50456118fb84da14366eceb94 Mon Sep 17 00:00:00 2001
From: Ziyue Zhang <Ziyue_Zhang@localhost.localdomain>
Date: Tue, 16 Mar 2021 15:41:34 -0500
Subject: [PATCH] WIP: Replacing old functions

---
 benchmarks/dvmp/analysis/dvmp.h       | 101 ++------------------------
 benchmarks/dvmp/analysis/vm_invar.cxx |  13 +---
 options/tracker_reconstruction.py     |   2 +-
 3 files changed, 13 insertions(+), 103 deletions(-)

diff --git a/benchmarks/dvmp/analysis/dvmp.h b/benchmarks/dvmp/analysis/dvmp.h
index 7271465d..f79ba0f2 100644
--- a/benchmarks/dvmp/analysis/dvmp.h
+++ b/benchmarks/dvmp/analysis/dvmp.h
@@ -19,28 +19,9 @@ namespace util {
 
   //========================================================================================================
   // for structure functions
-
   struct inv_quant { // add more when needed
     double nu, Q2, x, y, t;
   };
-
-  // for simu
-  /*inline inv_quant calc_inv_quant_sim(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]);//exact P' - P
-    //ROOT::Math::PxPyPzMVector Delta(parts[0] - parts[2] - parts[5]);//e - e' - jpsi
-    //ROOT::Math::PxPyPzMVector Delta(parts[0] - parts[2] - parts[7] - parts[8]);//jpsi->l l' gamma, ignore gamma
-    //P' - P = e - e' - jpsi, jpsi != l + l' (there could be an additional photon)
-    
-    
-    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;
-  }*/
   //import Geant4 and set the wanted particles in the intended order
   //0:e0  1:p0    2:e1    3:p1    4:recoil system (without p1)    5:l1 from 4 6:l2 from 4
   inline auto momenta_sort_sim(const std::vector<dd4pod::Geant4ParticleData>& parts, std::string_view mother, std::string_view daughter){//mother and daughter are not used yet; will be useful when generater is different and/or when the mcparticles doesn't follow the same order in all events
@@ -54,10 +35,10 @@ namespace util {
       double e = sqrt(px*px + py*py + pz*pz + mass*mass);
       momenta[i].SetPxPyPzE(px, py, pz, e);
     }
-    for(int i = 0 ; i < 7 ; i++){
-        cout<<Form("sim, idx = %d, P(px, py, pz, mass) = (%f, %f, %f, %f)", i, momenta[i].px(), momenta[i].py(), momenta[i].pz(), momenta[i].mass())<<endl;
-    }
-    cout<<"========================================="<<endl;
+    //for(int i = 0 ; i < 7 ; i++){
+    //    cout<<Form("sim, idx = %d, P(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;
   }
   
@@ -132,79 +113,13 @@ namespace util {
       //  dp = 10. - ptmp;
       //}
     }
-    for(int i = 0 ; i < 7 ; i++){
-        cout<<Form("dum, idx = %d, P(px, py, pz, mass) = (%f, %f, %f, %f)", i, momenta[i].px(), momenta[i].py(), momenta[i].pz(), momenta[i].mass())<<endl;
-    }
-    cout<<"========================================="<<endl;
+    //for(int i = 0 ; i < 7 ; i++){
+    //    cout<<Form("dum, idx = %d, P(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;
   }
   
- //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;
-    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) {
-      if( fabs(parts[i].mass() - daughter_mass)/daughter_mass > 0.01) continue;
-      for (int j = i + 1; j < parts.size(); ++j) {
-        if( fabs(parts[j].mass() - daughter_mass)/daughter_mass > 0.01) continue;
-        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 ){      //fewer than 3 candidates 
-      inv_quant quantities = {-99999., -99999., -99999., -99999.};
-      return quantities;
-    }
-    
-    //construct the beam kinematics
-    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"));
-    double e1_Energy = sqrt((1.305e-8 - 10.)*(1.305e-8 - 10.) + (0.0005109+9.888e-8)*(0.0005109+9.888e-8));
-    double P_Energy = sqrt((99.995598 + 1.313e-7)*(99.995598 + 1.313e-7) + (0.938272-1.23e-12)*(0.938272-1.23e-12));
-    e1.SetPxPyPzE(0., 0., 1.305e-8 - 10., e1_Energy);
-    P.SetPxPyPzE(0., 0., 99.995598 + 1.313e-7, P_Energy);
-    int scatteredIdx = -1;
-    
-    //float dp = 10.;
-    for(int i = 0 ; i < parts.size(); i++){
-      if(i==first || i==second) continue;   //skip the paired leptons
-      //if( fabs(parts[i].mass() - get_pdg_mass("electron"))/get_pdg_mass("electron") > 0.01) continue;
-      if( fabs(parts[i].mass() - 0.0005109 - 9.888e-8)/(0.0005109 + 9.888e-8) > 0.01) continue;
-      //ROOT::Math::PxPyPzMVector k_prime(parts[i]);  //scattered
-      //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 = {-99999., -99999., -99999., -99999.};
-      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};
-    return quantities;
-  }*/
-  
   inline inv_quant calc_inv_quant(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
diff --git a/benchmarks/dvmp/analysis/vm_invar.cxx b/benchmarks/dvmp/analysis/vm_invar.cxx
index 0a93c31a..df32f452 100644
--- a/benchmarks/dvmp/analysis/vm_invar.cxx
+++ b/benchmarks/dvmp/analysis/vm_invar.cxx
@@ -76,25 +76,20 @@ int vm_invar(const std::string& config_name)
   };*/
   auto momenta_sort_sim = [vm_name, decay_name](const std::vector<dd4pod::Geant4ParticleData>& parts){
     return util::momenta_sort_sim(parts, vm_name, decay_name);
-  };
   auto momenta_sort_rec = [vm_name, decay_name](const std::vector<eic::ReconstructedParticleData>& parts){
+  };
     return util::momenta_sort_rec(parts, vm_name, decay_name);
   };
   //====================================================================
 
   // Define analysis flow
-  auto d_im = d.Define("p_rec", util::momenta_RC, {"DummyReconstructedParticles"})
-                  .Define("p_sim", util::momenta_from_simulation, {"mcparticles2"})
+  auto d_im = d.//Define("p_rec", util::momenta_RC, {"DummyReconstructedParticles"})
+                  //.Define("p_sim", util::momenta_from_simulation, {"mcparticles2"})
                   .Define("p_rec_sorted", momenta_sort_rec, {"DummyReconstructedParticles"})
                   .Define("p_sim_sorted", momenta_sort_sim, {"mcparticles2"})
-                  .Define("N", "p_rec.size()")
-                  
-                  
+                  .Define("N", "p_rec_sorted.size()")
                   .Define("invariant_quantities_rec", util::calc_inv_quant, {"p_rec_sorted"})
                   .Define("invariant_quantities_sim", util::calc_inv_quant, {"p_sim_sorted"})
-                  
-                  
-                  
                   //.Define("p_sim", util::momenta_from_simulation, {"mcparticles2"})
                   //================================================================
                   //.Define("invariant_quantities_rec", calc_inv_quant_rec, {"p_rec"})
diff --git a/options/tracker_reconstruction.py b/options/tracker_reconstruction.py
index 6bef5d7b..e303b38e 100644
--- a/options/tracker_reconstruction.py
+++ b/options/tracker_reconstruction.py
@@ -57,7 +57,7 @@ podioinput = PodioInput("PodioReader",
 dummy = MC2DummyParticle("MC2Dummy",
         inputCollection="mcparticles",
         outputCollection="DummyReconstructedParticles",
-        smearing = 0.0)
+        smearing = 0.1)
 
 ## copiers to get around input --> output copy bug. Note the "2" appended to the output collection.
 copier = MCCopier("MCCopier", 
-- 
GitLab