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