File indexing completed on 2024-09-27 07:02:59
0001
0002
0003
0004 #include <edm4eic/EDM4eicVersion.h>
0005 #if EDM4EIC_VERSION_MAJOR >= 6
0006
0007 #include <Math/GenVector/LorentzVector.h>
0008 #include <Math/GenVector/PxPyPzE4D.h>
0009 #include <Math/Vector4Dfwd.h>
0010 #include <edm4eic/InclusiveKinematicsCollection.h>
0011 #include <edm4eic/ReconstructedParticleCollection.h>
0012 #include <edm4hep/MCParticleCollection.h>
0013 #include <edm4hep/Vector3f.h>
0014 #include <fmt/core.h>
0015 #include <cmath>
0016 #include <gsl/pointers>
0017
0018 #include "Beam.h"
0019 #include "Boost.h"
0020 #include "InclusiveKinematicsESigma.h"
0021
0022 using ROOT::Math::PxPyPzEVector;
0023
0024 namespace eicrecon {
0025
0026 void InclusiveKinematicsESigma::init() { }
0027
0028 void InclusiveKinematicsESigma::process(
0029 const InclusiveKinematicsESigma::Input& input,
0030 const InclusiveKinematicsESigma::Output& output) const {
0031
0032 const auto [mcparts, escat, hfs] = input;
0033 auto [kinematics] = output;
0034
0035
0036 const auto ei_coll = find_first_beam_electron(mcparts);
0037 if (ei_coll.size() == 0) {
0038 debug("No beam electron found");
0039 return;
0040 }
0041 const PxPyPzEVector ei(
0042 round_beam_four_momentum(
0043 ei_coll[0].getMomentum(),
0044 m_particleSvc.particle(ei_coll[0].getPDG()).mass,
0045 {-5.0, -10.0, -18.0},
0046 0.0)
0047 );
0048
0049
0050 const auto pi_coll = find_first_beam_hadron(mcparts);
0051 if (pi_coll.size() == 0) {
0052 debug("No beam hadron found");
0053 return;
0054 }
0055 const PxPyPzEVector pi(
0056 round_beam_four_momentum(
0057 pi_coll[0].getMomentum(),
0058 m_particleSvc.particle(pi_coll[0].getPDG()).mass,
0059 {41.0, 100.0, 275.0},
0060 m_crossingAngle)
0061 );
0062
0063
0064 auto boost = determine_boost(ei, pi);
0065
0066
0067 auto kf = escat->at(0);
0068 PxPyPzEVector e_lab(kf.getMomentum().x, kf.getMomentum().y, kf.getMomentum().z, kf.getEnergy());
0069 PxPyPzEVector e_boosted = apply_boost(boost, e_lab);
0070 auto pt_e = e_boosted.Pt();
0071 auto sigma_e = e_boosted.E() - e_boosted.Pz();
0072
0073
0074 auto sigma_h = hfs->at(0).getSigma();
0075 auto ptsum = hfs->at(0).getPT();
0076 auto gamma_h = hfs->at(0).getGamma();
0077
0078 if (sigma_h <= 0) {
0079 debug("No scattered electron found or sigma zero or negative");
0080 return;
0081 }
0082
0083 auto sigma_tot = sigma_e + sigma_h;
0084
0085
0086 const auto y_e = 1. - sigma_e / (2.*ei.energy());
0087 const auto Q2_e = (pt_e*pt_e) / (1. - y_e);
0088
0089 const auto y_sig = sigma_h / sigma_tot;
0090 const auto Q2_sig = (pt_e*pt_e) / (1. - y_sig);
0091 const auto x_sig = Q2_sig / (4.*ei.energy()*pi.energy()*y_sig);
0092
0093 static const auto m_proton = m_particleSvc.particle(2212).mass;
0094 const auto Q2_esig = Q2_e;
0095 const auto x_esig = x_sig;
0096 const auto y_esig = Q2_esig / (4.*ei.energy()*pi.energy()*x_esig);
0097 const auto nu_esig = Q2_esig / (2.*m_proton*x_esig);
0098 const auto W_esig = sqrt(m_proton*m_proton + 2*m_proton*nu_esig - Q2_esig);
0099 auto kin = kinematics->create(x_esig, Q2_esig, W_esig, y_esig, nu_esig);
0100 kin.setScat(kf);
0101
0102
0103 debug("x,Q2,W,y,nu = {},{},{},{},{}", kin.getX(),
0104 kin.getQ2(), kin.getW(), kin.getY(), kin.getNu());
0105 }
0106
0107 }
0108 #endif