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