change

benchmarks/dvcs-d/analysis/dvcs_d_analysis.cxx

@@ -69,6 +69,7 @@ int dvcs_d_analysis(const std::string& config_name)
              .Define("gammaREC",findGamma,{"ReconstructedParticles"}).Define("MassREC",getMass,{"gammaREC"})
              .Define("gamma_rec_pt",getPt,{"gammaREC"}).Define("gamma_rec_eta",getEta,{"gammaREC"}).Define("gamma_rec_E",getE,{"gammaREC"})
              .Define("protonREC",findScatProton,{"ReconstructedFFParticles"}).Define("proton_rec_eta",getEta,{"protonREC"}).Define("proton_rec_phi",getPhi,{"protonREC"}).Define("proton_rec_px",getPx,{"protonREC"}).Define("proton_rec_py",getPy,{"protonREC"}).Define("proton_rec_momentum",getPmiss,{"protonREC","mcparticles"})
+             .Define("proton_rec_p",getP,{"protonREC"}).Define("proton_rec_pz",getPz,{"protonREC"})
              .Define("neutronREC",findScatNeutron,{"ReconstructedFFParticles"}).Define("neutron_rec_px",getPx,{"neutronREC"}).Define("neutron_rec_py",getPy,{"neutronREC"})
              .Define("proton_rec_pt",getPt,{"protonREC"}).Define("proton_rec_theta",getTheta,{"protonREC"})
              .Define("t_REC",giveme_t_doubleTagging_REC,{"ReconstructedFFParticles","mcparticles"}).Define("t_REC_A",giveme_t,{"gammaREC","scatElec"})
@@ -85,6 +86,8 @@ int dvcs_d_analysis(const std::string& config_name)
   auto h_Eta_proton_REC = d1.Histo1D({"h_Eta_proton_REC", "; #eta; counts", 100, 0,9}, "proton_rec_eta");
   auto h_Theta_proton_REC = d1.Histo1D({"h_Theta_proton_REC", "; #theta; counts", 100, 0,0.1}, "proton_rec_theta");
   auto h_Momentum_proton_REC = d1.Histo1D({"h_Momentum_proton_REC", "; P; counts", 100, 0, 2}, "proton_rec_momentum");
+  auto h_P_proton_REC = d1.Histo1D({"h_P_proton_REC", "; P; counts", 200, 0, 200}, "proton_rec_p");
+  auto h_Pz_proton_REC = d1.Histo1D({"h_Pz_proton_REC", "; P; counts", 200, 0, 200}, "proton_rec_pz");
   auto h_EtaPhi_proton_REC = d1.Histo2D({"h_EtaPhi_proton_REC",";#eta;#phi",100,-10,10,100,-PI,PI},"proton_rec_eta","proton_rec_phi");
   auto h_PxPy_proton_REC = d1.Histo2D({"h_PxPy_proton_REC",";#px;#py",200,-6,2,200,-6,2},"proton_rec_px","proton_rec_py");
   auto h_PxPy_neutron_REC = d1.Histo2D({"h_PxPy_neutron_REC",";#px;#py",200,-6,2,200,-6,2},"neutron_rec_px","neutron_rec_py");
@@ -102,7 +105,7 @@ int dvcs_d_analysis(const std::string& config_name)
              .Define("scatElecMC",findScatElecMC, {"mcparticles"}).Define("etaElecMC",getEta,{"scatElecMC"})
              .Define("gammaMC",findGammaMC,{"mcparticles"}).Define("MassMC",getMass,{"gammaMC"}).Define("gamma_mc_pt",getPt,{"gammaMC"}).Define("gamma_mc_eta",getEta,{"gammaMC"})
              .Define("protonMC",findScatProtonMC,{"mcparticles"}).Define("proton_mc_pt",getPt,{"protonMC"}).Define("proton_mc_eta",getEta,{"protonMC"}).Define("proton_mc_theta",getTheta,{"protonMC"})
-             .Define("proton_mc_momentum",getPmiss,{"protonMC","mcparticles"})
+             .Define("proton_mc_momentum",getPmiss,{"protonMC","mcparticles"}).Define("proton_mc_p",getP,{"protonMC"}).Define("proton_mc_pz",getPz,{"protonMC"})
              .Define("t_MC",giveme_t_neutron_MC,{"mcparticles"})
              .Filter(kineCut,{"Q2_elec","y_elec"});
 
@@ -114,6 +117,8 @@ int dvcs_d_analysis(const std::string& config_name)
   auto h_Eta_proton_MC = d2.Histo1D({"h_Eta_proton_MC", "; #eta; counts", 100, 0, 9}, "proton_mc_eta");
   auto h_Theta_proton_MC = d2.Histo1D({"h_Theta_proton_MC", "; #theta; counts", 100, 0, 0.1}, "proton_mc_theta");
   auto h_Momentum_proton_MC = d2.Histo1D({"h_Momentum_proton_MC", "; P; counts", 100, 0, 2}, "proton_mc_momentum");
+  auto h_P_proton_MC = d2.Histo1D({"h_P_proton_MC", "; P; counts", 200, 0, 200}, "proton_mc_p");
+  auto h_Pz_proton_MC = d2.Histo1D({"h_Pz_proton_MC", "; P; counts", 200, 0, 200}, "proton_mc_pz");
   auto h_t_MC = d2.Histo1D({"h_t_MC", "; ; counts", 50, 0, 2}, "t_MC");
  
   TString output_name_dir = output_prefix.c_str();
@@ -141,6 +146,8 @@ int dvcs_d_analysis(const std::string& config_name)
   h_Eta_proton_REC->Write();
   h_Theta_proton_REC->Write();
   h_Momentum_proton_REC->Write();
+  h_P_proton_REC->Write();
+  h_Pz_proton_REC->Write();
   h_EtaPhi_proton_REC->Write();
   h_PxPy_proton_REC->Write();
   h_PxPy_neutron_REC->Write();
@@ -157,6 +164,8 @@ int dvcs_d_analysis(const std::string& config_name)
   h_Eta_proton_MC->Write();
   h_Theta_proton_MC->Write();
   h_Momentum_proton_MC->Write();
+  h_P_proton_MC->Write();
+  h_Pz_proton_MC->Write();
   h_t_MC->Write();
 
   output->Write();

benchmarks/dvcs-d/analysis/pleaseIncludeMe.h

@@ -330,6 +330,30 @@ auto getPmiss(const std::vector<ROOT::Math::PxPyPzMVector>& mom,
   return momentumVec;
 }
 
+auto getP(const std::vector<ROOT::Math::PxPyPzMVector>& mom)
+{
+  std::vector<double> pVec;
+  for(auto& i1:mom){
+    double momentum = i1.P();
+    if(i1.Px()<-1e9){momentum=-10.;}
+    pVec.push_back(momentum);
+  }
+  return pVec;
+
+}
+
+auto getPz(const std::vector<ROOT::Math::PxPyPzMVector>& mom)
+{
+  std::vector<double> pzVec;
+  for(auto& i1:mom){
+    double pz = i1.Pz();
+    if(i1.Px()<-1e9){pz=-10.;}
+    pzVec.push_back(pz);
+  }
+  return pzVec;
+
+}
+
 auto getAngleDiff(const std::vector<ROOT::Math::PxPyPzMVector> ph_gen,
   const std::vector<ROOT::Math::PxPyPzMVector> ph_rec)
 {

options/reconstruction.py

@@ -647,7 +647,7 @@ fast_ff = FFSmearedParticles("fast_ff",
         enableB0         = True,
         enableRP         = True,
         enableOMD        = True,
-        ionBeamEnergy    = 200,
+        ionBeamEnergy    = 100,
         crossingAngle    = -0.025)
 algorithms.append(fast_ff)