algorithms/reco/InclusiveKinematicsSigma.cc

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 Wouter Deconinck, Barak Schmookler
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 <edm4hep/MCParticleCollection.h>
0012 #include <edm4hep/Vector3f.h>
0013 #include <fmt/core.h>
0014 #include <cmath>
0015 #include <gsl/pointers>
0016
0017 #include "Beam.h"
0018 #include "Boost.h"
0019 #include "InclusiveKinematicsSigma.h"
0020
0021 using ROOT::Math::PxPyPzEVector;
0022
0023 namespace eicrecon {
0024
0025   void InclusiveKinematicsSigma::init() { }
0026
0027   void InclusiveKinematicsSigma::process(
0028       const InclusiveKinematicsSigma::Input& input,
0029       const InclusiveKinematicsSigma::Output& output) const {
0030
0031     const auto [mcparts, escat, hfs] = input;
0032     auto [kinematics] = output;
0033
0034     // Get incoming electron beam
0035     const auto ei_coll = find_first_beam_electron(mcparts);
0036     if (ei_coll.size() == 0) {
0037       debug("No beam electron found");
0038       return;
0039     }
0040     const PxPyPzEVector ei(
0041       round_beam_four_momentum(
0042         ei_coll[0].getMomentum(),
0043         m_particleSvc.particle(ei_coll[0].getPDG()).mass,
0044         {-5.0, -10.0, -18.0},
0045         0.0)
0046       );
0047
0048     // Get incoming hadron beam
0049     const auto pi_coll = find_first_beam_hadron(mcparts);
0050     if (pi_coll.size() == 0) {
0051       debug("No beam hadron found");
0052       return;
0053     }
0054     const PxPyPzEVector pi(
0055       round_beam_four_momentum(
0056         pi_coll[0].getMomentum(),
0057         m_particleSvc.particle(pi_coll[0].getPDG()).mass,
0058         {41.0, 100.0, 275.0},
0059         m_crossingAngle)
0060       );
0061
0062     // Get boost to colinear frame
0063     auto boost = determine_boost(ei, pi);
0064
0065     // Get electron variables
0066     auto kf = escat->at(0);
0067     PxPyPzEVector e_lab(kf.getMomentum().x, kf.getMomentum().y, kf.getMomentum().z, kf.getEnergy());
0068     PxPyPzEVector e_boosted = apply_boost(boost, e_lab);
0069     auto pt_e = e_boosted.Pt();
0070     auto sigma_e = e_boosted.E() - e_boosted.Pz();
0071
0072     // Get hadronic final state variables
0073     auto sigma_h = hfs->at(0).getSigma();
0074     auto ptsum = hfs->at(0).getPT();
0075     auto gamma_h = hfs->at(0).getGamma();
0076
0077     if (sigma_h <= 0) {
0078       debug("No scattered electron found or sigma zero or negative");
0079       return;
0080     }
0081
0082     auto sigma_tot = sigma_e + sigma_h;
0083
0084     // Calculate kinematic variables
0085     static const auto m_proton = m_particleSvc.particle(2212).mass;
0086     const auto y_sig = sigma_h / sigma_tot;
0087     const auto Q2_sig = (pt_e*pt_e) / (1. - y_sig);
0088     const auto x_sig = Q2_sig / (4.*ei.energy()*pi.energy()*y_sig);
0089     const auto nu_sig = Q2_sig / (2.*m_proton*x_sig);
0090     const auto W_sig = sqrt(m_proton*m_proton + 2*m_proton*nu_sig - Q2_sig);
0091     auto kin = kinematics->create(x_sig, Q2_sig, W_sig, y_sig, nu_sig);
0092     kin.setScat(kf);
0093
0094     debug("x,Q2,W,y,nu = {},{},{},{},{}", kin.getX(),
0095             kin.getQ2(), kin.getW(), kin.getY(), kin.getNu());
0096   }
0097
0098 }
0099 #endif