change dvcs-d analysis

benchmarks/dvcs-d/analysis/dvcs_d_analysis.cxx

@@ -70,7 +70,7 @@ int dvcs_d_analysis(const std::string& config_name)
              .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_pt",getPt,{"protonREC"}).Define("proton_rec_theta",getTheta,{"protonREC"})
-             .Define("t_REC",giveme_t_REC,{"protonREC","mcparticles"}).Define("t_REC_A",giveme_t,{"gammaREC","scatElec"})
+             .Define("t_REC",giveme_t_doubleTagging_REC,{"ReconstructedFFParticles"}).Define("t_REC_A",giveme_t,{"gammaREC","scatElec"})
              .Filter(kineCut,{"Q2_elec","y_elec"});
 
   auto h_Q2_elec = d1.Histo1D({"h_Q2_elec", "; GeV^2; counts", 100, -5, 25}, "Q2_elec");
@@ -98,7 +98,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("t_MC",giveme_t_MC,{"mcparticles"})
+             .Define("t_MC",giveme_t_neutron_MC,{"mcparticles"})
              .Filter(kineCut,{"Q2_elec","y_elec"});
 
   auto h_Eta_scatElec_MC = d2.Histo1D({"h_Eta_scatElec_MC",";eta; counts",100,-11,9}, "etaElecMC");

benchmarks/dvcs-d/analysis/pleaseIncludeMe.h

@@ -311,6 +311,24 @@ auto giveme_t_MC(const std::vector<dd4pod::Geant4ParticleData>& parts){
   return t_vec;
 }
 
+auto giveme_t_neutron_MC(const std::vector<dd4pod::Geant4ParticleData>& parts){
+  std::vector<double > t_vec;
+  TLorentzVector pIn,pOut;
+  for(auto& i1 : parts){
+    if(i1.genStatus==4&&i1.pdgID==2112) {
+      TVector3 pIn_v3(i1.ps.x,i1.ps.y,i1.ps.z);
+      pIn.SetVectM(pIn_v3,i1.mass);
+    }
+    if(i1.genStatus==1&&i1.pdgID==2112) {
+      TVector3 pOut_v3(i1.ps.x,i1.ps.y,i1.ps.z);
+      pOut.SetVectM(pOut_v3,i1.mass);
+    }
+  }
+  t_vec.push_back( -(pOut-pIn).Mag2() );
+  
+  return t_vec;
+}
+
 auto giveme_t_REC(const std::vector<ROOT::Math::PxPyPzMVector>& mom,
   const std::vector<dd4pod::Geant4ParticleData>& parts){
   
@@ -331,6 +349,25 @@ auto giveme_t_REC(const std::vector<ROOT::Math::PxPyPzMVector>& mom,
   return t_vec;
 }
 
+auto giveme_t_doubleTagging_REC(const std::vector<ROOT::Math::PxPyPzMVector>& mom){
+  
+  std::vector<double> t_vec;
+  TLorentzVector nOut,pOut;
+  for(auto&i2: mom){
+    if(fabs(i2.mass-0.93957)<1e-4){
+      TVector3 nOut_v3(i2.Px(),i2.Py(),i2.Pz());
+      nOut.SetVectM(nOut_v3,0.93957);
+    }
+    if(fabs(i2.mass-0.93827)<1e-4){
+      TVector3 pOut_v3(-i2.Px(),-i2.Py(),-i2.Pz());//remember the - sign.
+      pOut.SetVectM(pOut_v3,0.93827);
+    }
+  }
+  t_vec.push_back( -(nOut-pOut).Mag2() );
+  
+  return t_vec;
+}
+
 auto giveme_t = [](std::vector<ROOT::Math::PxPyPzMVector> vm, 
    std::vector<ROOT::Math::PxPyPzMVector> scatElec){
   std::vector<double > t_vec;