algorithms/reco/TransformBreitFrame.cc

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2024 John Lajoie
0003
0004 // class definition
0005 #include "TransformBreitFrame.h"
0006
0007 #include <Math/GenVector/Boost.h>
0008 #include <Math/GenVector/Cartesian3D.h>
0009 #include <Math/GenVector/DisplacementVector3D.h>
0010 #include <Math/GenVector/LorentzVector.h>
0011 #include <Math/GenVector/PxPyPzE4D.h>
0012 #include <Math/GenVector/Rotation3D.h>
0013 #include <Math/Vector4Dfwd.h>
0014 #include <edm4eic/Cov4f.h>
0015 #include <edm4eic/Vertex.h>
0016 #include <edm4hep/Vector3f.h>
0017 #include <edm4hep/utils/kinematics.h>
0018 #include <fmt/core.h>
0019 #include <gsl/pointers>
0020 #include <vector>
0021
0022 #include "Beam.h"
0023
0024 namespace eicrecon {
0025
0026 void TransformBreitFrame::process(const TransformBreitFrame::Input& input,
0027                                   const TransformBreitFrame::Output& output) const {
0028   // Grab input collections
0029   const auto [mcpart, kine, lab_collection] = input;
0030   auto [breit_collection]                   = output;
0031
0032   // Beam momenta extracted from MCParticle
0033   // This is the only place truth information is used!
0034
0035   // Get incoming electron beam
0036   const auto ei_coll = find_first_beam_electron(mcpart);
0037   if (ei_coll.empty()) {
0038     debug("No beam electron found");
0039     return;
0040   }
0041   const PxPyPzEVector e_initial(round_beam_four_momentum(
0042       ei_coll[0].getMomentum(), m_particleSvc.particle(ei_coll[0].getPDG()).mass,
0043       {-5.0, -10.0, -18.0}, 0.0));
0044
0045   // Get incoming hadron beam
0046   const auto pi_coll = find_first_beam_hadron(mcpart);
0047   if (pi_coll.empty()) {
0048     debug("No beam hadron found");
0049     return;
0050   }
0051   const PxPyPzEVector p_initial(round_beam_four_momentum(
0052       pi_coll[0].getMomentum(), m_particleSvc.particle(pi_coll[0].getPDG()).mass,
0053       {41.0, 100.0, 275.0}, m_crossingAngle));
0054
0055   debug("electron energy, proton energy = {},{}", e_initial.E(), p_initial.E());
0056
0057   // Get the event kinematics, set up transform
0058   if (kine->empty()) {
0059     debug("No kinematics found");
0060     return;
0061   }
0062
0063   const auto& evt_kin = kine->at(0);
0064
0065   const auto meas_x  = evt_kin.getX();
0066   const auto meas_Q2 = evt_kin.getQ2();
0067
0068   // Use relation to get reconstructed scattered electron
0069   const PxPyPzEVector e_final =
0070       edm4hep::utils::detail::p4(evt_kin.getScat(), &edm4hep::utils::UseEnergy);
0071   debug("scattered electron in lab frame px,py,pz,E = {},{},{},{}", e_final.Px(), e_final.Py(),
0072         e_final.Pz(), e_final.E());
0073
0074   // Set up the transformation
0075   const PxPyPzEVector virtual_photon = (e_initial - e_final);
0076   debug("virtual photon in lab frame px,py,pz,E = {},{},{},{}", virtual_photon.Px(),
0077         virtual_photon.Py(), virtual_photon.Pz(), virtual_photon.E());
0078
0079   debug("x, Q^2 = {},{}", meas_x, meas_Q2);
0080
0081   // Set up the transformation (boost) to the Breit frame
0082   const auto P3 = p_initial.Vect();
0083   const auto q3 = virtual_photon.Vect();
0084   const ROOT::Math::Boost breit(-(2.0 * meas_x * P3 + q3) *
0085                                 (1.0 / (2.0 * meas_x * p_initial.E() + virtual_photon.E())));
0086
0087   PxPyPzEVector p_initial_breit      = (breit * p_initial);
0088   PxPyPzEVector e_initial_breit      = (breit * e_initial);
0089   PxPyPzEVector e_final_breit        = (breit * e_final);
0090   PxPyPzEVector virtual_photon_breit = (breit * virtual_photon);
0091
0092   // Now rotate so the virtual photon momentum is all along the negative z-axis
0093   const auto zhat = -virtual_photon_breit.Vect().Unit();
0094   const auto yhat = (e_initial_breit.Vect().Cross(e_final_breit.Vect())).Unit();
0095   const auto xhat = yhat.Cross(zhat);
0096
0097   const ROOT::Math::Rotation3D breitRotInv(xhat, yhat, zhat);
0098   const ROOT::Math::Rotation3D breitRot = breitRotInv.Inverse();
0099
0100   // Diagnostics
0101   p_initial_breit      = breitRot * p_initial_breit;
0102   e_initial_breit      = breitRot * e_initial_breit;
0103   e_final_breit        = breitRot * e_final_breit;
0104   virtual_photon_breit = breitRot * virtual_photon_breit;
0105
0106   debug("incoming hadron in Breit frame px,py,pz,E = {},{},{},{}", p_initial_breit.Px(),
0107         p_initial_breit.Py(), p_initial_breit.Pz(), p_initial_breit.E());
0108   debug("virtual photon in Breit frame px,py,pz,E = {},{},{},{}", virtual_photon_breit.Px(),
0109         virtual_photon_breit.Py(), virtual_photon_breit.Pz(), virtual_photon_breit.E());
0110
0111   // look over the input particles and transform
0112   for (const auto& lab : *lab_collection) {
0113
0114     // Transform to Breit frame
0115     PxPyPzEVector lab_particle(lab.getMomentum().x, lab.getMomentum().y, lab.getMomentum().z,
0116                                lab.getEnergy());
0117     PxPyPzEVector breit_particle = breitRot * (breit * lab_particle);
0118
0119     // create particle to store in output collection
0120     auto breit_out = breit_collection->create();
0121     breit_out.setMomentum(
0122         edm4hep::Vector3f(breit_particle.Px(), breit_particle.Py(), breit_particle.Pz()));
0123     breit_out.setEnergy(breit_particle.E());
0124
0125     // Copy the rest of the particle information
0126     breit_out.setType(lab.getType());
0127     breit_out.setReferencePoint(lab.getReferencePoint());
0128     breit_out.setCharge(lab.getCharge());
0129     breit_out.setMass(lab.getMass());
0130     breit_out.setGoodnessOfPID(lab.getGoodnessOfPID());
0131     breit_out.setCovMatrix(lab.getCovMatrix());
0132     breit_out.setPDG(lab.getPDG());
0133     breit_out.setStartVertex(lab.getStartVertex());
0134     breit_out.setParticleIDUsed(lab.getParticleIDUsed());
0135
0136     // set up a relation between the lab and Breit frame representations
0137     breit_out.addToParticles(lab);
0138   }
0139
0140 } // end 'process'
0141
0142 } // end namespace eicrecon