ElectronID/Boost/Boost.h

0001 // Modified from EICrecon Boost.h
0002
0003 // SPDX-License-Identifier: LGPL-3.0-or-later
0004 // Copyright (C) 2022 Wouter Deconinck, Barak Schmookler
0005
0006 #pragma once
0007
0008 #include <Math/Vector4D.h>
0009 #include <Math/LorentzRotation.h>
0010 #include <Math/LorentzVector.h>
0011 #include <Math/RotationX.h>
0012 #include <Math/RotationY.h>
0013 #include <Math/Boost.h>
0014
0015 using ROOT::Math::PxPyPzEVector;
0016
0017 namespace eicrecon {
0018
0019 using ROOT::Math::LorentzRotation;
0020
0021 inline LorentzRotation determine_boost(PxPyPzEVector ei, PxPyPzEVector pi) {
0022
0023   using ROOT::Math::Boost;
0024   using ROOT::Math::RotationX;
0025   using ROOT::Math::RotationY;
0026
0027   // Step 1: Find the needed boosts and rotations from the incoming lepton and hadron beams
0028   // (note, this will give you a perfect boost, in principle you will not know the beam momenta exactly and should use an average)
0029
0030   PxPyPzEVector eo = ei;
0031   PxPyPzEVector po = pi;
0032
0033   // cout << "lab e and p: " << Form("(%f, %f, %f, %f) (%f, %f, %f, %f)", ei.Px(), ei.Py(), ei.Pz(), ei.E(), pi.Px(), pi.Py(), pi.Pz(), pi.E()) << endl;
0034
0035   // Define the Boost to make beams back-to-back
0036   const auto cmBoost = (ei + pi).BoostToCM();
0037   const Boost boost_to_cm(cmBoost);
0038
0039   // Boost to COM frame
0040   pi = boost_to_cm(pi);
0041   ei = boost_to_cm(ei);
0042   // cout << "CM boosted e and p: " << Form("(%f, %f, %f, %f) (%f, %f, %f, %f)", ei.Px(), ei.Py(), ei.Pz(), ei.E(), pi.Px(), pi.Py(), pi.Pz(), pi.E()) << endl;
0043
0044   // This will boost beams from a center of momentum frame back to (nearly) their original energies
0045   PxPyPzEVector eh(0, 0, -1*sqrt(pow(eo.E(),2)-pow(eo.M(),2)), eo.E());
0046   PxPyPzEVector ph(0, 0,    sqrt(pow(po.E(),2)-pow(po.M(),2)), po.E());
0047   // cout << "headon e and p: " << Form("(%f, %f, %f, %f) (%f, %f, %f, %f)", eh.Px(), eh.Py(), eh.Pz(), eh.E(), ph.Px(), ph.Py(), ph.Pz(), ph.E()) << endl;
0048
0049   const auto hoBoost = (eh + ph).BoostToCM();
0050   // const Boost headon_to_cm(hoBoost);
0051   const Boost boost_to_headon(-hoBoost);
0052
0053   // PxPyPzEVector et = headon_to_cm(eh);
0054   // PxPyPzEVector pt = headon_to_cm(ph);
0055   // cout << "headon boosted to CM e and p: " << Form("(%f, %f, %f, %f) (%f, %f, %f, %f)", et.Px(), et.Py(), et.Pz(), et.E(), pt.Px(), pt.Py(), pt.Pz(), pt.E()) << endl;
0056
0057   // PxPyPzEVector erb = boost_to_headon(et);
0058   // PxPyPzEVector prb = boost_to_headon(pt);
0059   // cout << "reversed boost e and p: " << Form("(%f, %f, %f, %f) (%f, %f, %f, %f)", erb.Px(), erb.Py(), erb.Pz(), erb.E(), prb.Px(), prb.Py(), prb.Pz(), prb.E()) << endl;
0060
0061   // Boost and rotate the incoming beams to find the proper rotations TLorentzVector
0062
0063   // Rotate to head-on
0064   RotationY rotAboutY(-1.0 * atan2(pi.Px(), pi.Pz())); // Rotate to remove x component of beams
0065   RotationX rotAboutX(+1.0 * atan2(pi.Py(), pi.Pz())); // Rotate to remove y component of beams
0066
0067   // PxPyPzEVector er = rotAboutX(rotAboutY(ei));
0068   // PxPyPzEVector pr = rotAboutX(rotAboutY(pi));
0069   // cout << "rotation e and p: " << Form("(%f, %f, %f, %f) (%f, %f, %f, %f)", er.Px(), er.Py(), er.Pz(), er.E(), pr.Px(), pr.Py(), pr.Pz(), pr.E()) << endl;
0070
0071   // PxPyPzEVector ef = boost_to_headon(er);
0072   // PxPyPzEVector pf = boost_to_headon(pr);
0073   // cout << "final boost e and p: " << Form("(%f, %f, %f, %f) (%f, %f, %f, %f)", ef.Px(), ef.Py(), ef.Pz(), ef.E(), pf.Px(), pf.Py(), pf.Pz(), pf.E()) << endl;
0074
0075   // cout << "**" << endl;
0076
0077   // final matrix: P' = [BtoH][RX][RY][BtoCM]P <-- Matrix multi. goes from R to L
0078   LorentzRotation tf;
0079   tf *= boost_to_headon;
0080   tf *= rotAboutX;
0081   tf *= rotAboutY;
0082   tf *= boost_to_cm;
0083
0084   // PxPyPzEVector em = tf(eo);
0085   // PxPyPzEVector pm = tf(po);
0086   // cout << "boost matrix e and p: " << Form("(%f, %f, %f, %f) (%f, %f, %f, %f)", em.Px(), em.Py(), em.Pz(), em.E(), pm.Px(), pm.Py(), pm.Pz(), pm.E()) << endl;
0087
0088   return tf;
0089 }
0090
0091 inline PxPyPzEVector apply_boost(const LorentzRotation& tf, PxPyPzEVector part) {
0092
0093   // Step 2: Apply boosts and rotations to any particle 4-vector
0094   // (here too, choices will have to be made as to what the 4-vector is for reconstructed particles)
0095
0096   // Boost and rotate particle 4-momenta into the headon frame
0097   tf(part);
0098   return part;
0099 }
0100
0101 } // namespace eicrecon