algorithms/reco/InclusiveKinematicsDA.cc

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 Wouter Deconinck
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 "InclusiveKinematicsDA.h"
0021
0022 using ROOT::Math::PxPyPzEVector;
0023
0024 namespace eicrecon {
0025
0026   void InclusiveKinematicsDA::init() {
0027     // m_pidSvc = service("ParticleSvc");
0028     // if (!m_pidSvc) {
0029     //   m_log->debug("Unable to locate Particle Service. "
0030     //     "Make sure you have ParticleSvc in the configuration.");
0031     // }
0032   }
0033
0034   void InclusiveKinematicsDA::process(
0035       const InclusiveKinematicsDA::Input& input,
0036       const InclusiveKinematicsDA::Output& output) const {
0037
0038     const auto [mcparts, escat, hfs] = input;
0039     auto [kinematics] = output;
0040
0041     // Get incoming electron beam
0042     const auto ei_coll = find_first_beam_electron(mcparts);
0043     if (ei_coll.size() == 0) {
0044       debug("No beam electron found");
0045       return;
0046     }
0047     const PxPyPzEVector ei(
0048       round_beam_four_momentum(
0049         ei_coll[0].getMomentum(),
0050         m_electron,
0051         {-5.0, -10.0, -18.0},
0052         0.0)
0053       );
0054
0055     // Get incoming hadron beam
0056     const auto pi_coll = find_first_beam_hadron(mcparts);
0057     if (pi_coll.size() == 0) {
0058       debug("No beam hadron found");
0059       return;
0060     }
0061     const PxPyPzEVector pi(
0062       round_beam_four_momentum(
0063         pi_coll[0].getMomentum(),
0064         pi_coll[0].getPDG() == 2212 ? m_proton : m_neutron,
0065         {41.0, 100.0, 275.0},
0066         m_crossingAngle)
0067       );
0068
0069     // Get boost to colinear frame
0070     auto boost = determine_boost(ei, pi);
0071
0072     // Get electron angle
0073     auto kf = escat->at(0);
0074     PxPyPzEVector e_lab(kf.getMomentum().x, kf.getMomentum().y, kf.getMomentum().z, kf.getEnergy());
0075     PxPyPzEVector e_boosted = apply_boost(boost, e_lab);
0076     auto theta_e = e_boosted.Theta();
0077
0078     // Get hadronic final state variables
0079     auto sigma_h = hfs->at(0).getSigma();
0080     auto ptsum = hfs->at(0).getPT();
0081     auto gamma_h = hfs->at(0).getGamma();
0082
0083     // Sigma zero or negative
0084     if (sigma_h <= 0) {
0085       debug("Sigma zero or negative");
0086       return;
0087     }
0088
0089     // Calculate kinematic variables
0090     const auto y_da = tan(gamma_h/2.) / ( tan(theta_e/2.) + tan(gamma_h/2.) );
0091     const auto Q2_da = 4.*ei.energy()*ei.energy() * ( 1. / tan(theta_e/2.) ) * ( 1. / (tan(theta_e/2.) + tan(gamma_h/2.)) );
0092     const auto x_da = Q2_da / (4.*ei.energy()*pi.energy()*y_da);
0093     const auto nu_da = Q2_da / (2.*m_proton*x_da);
0094     const auto W_da = sqrt(m_proton*m_proton + 2*m_proton*nu_da - Q2_da);
0095     auto kin = kinematics->create(x_da, Q2_da, W_da, y_da, nu_da);
0096     kin.setScat(kf);
0097
0098     debug("x,Q2,W,y,nu = {},{},{},{},{}", kin.getX(),
0099             kin.getQ2(), kin.getW(), kin.getY(), kin.getNu());
0100   }
0101
0102 }
0103 #endif