include/jpacPhoto/reaction_kinematics.hpp

0001 // Class to contain all relevant kinematic quantities. The kinematics of the reaction

0002 // gamma p -> X p' is entirely determined by specifying the mass of the vector particle.

0003 //

0004 // Additional options to include virtual photon and different baryons (e.g. gamma p -> X Lambda_c)

0005 // also available

0006 //

0007 // Author:       Daniel Winney (2020)

0008 // Affiliation:  Joint Physics Analysis Center (JPAC)

0009 // Email:        dwinney@iu.edu

0010 // ---------------------------------------------------------------------------

0011
0012
0013 #ifndef _KINEMATICS_
0014 #define _KINEMATICS_
0015
0016 #include "constants.hpp"
0017 #include "misc_math.hpp"
0018 #include "two_body_state.hpp"
0019 #include "dirac_spinor.hpp"
0020 #include "polarization_vector.hpp"
0021 #include "helicities.hpp"
0022
0023 #include "TMath.h"
0024
0025 #include <array>
0026 #include <vector>
0027 #include <string>
0028 #include <cmath>
0029
0030 namespace jpacPhoto
0031 {
0032     // ---------------------------------------------------------------------------

0033     // The reaction kinematics object is intended to have all relevant kinematic quantities

0034     // forthe reaction. Here you'll find the momenta and energies of all particles,

0035     //  spinors for the baryons and polarization vectors for the gamma and produced meson

0036     // ---------------------------------------------------------------------------

0037
0038     class reaction_kinematics
0039     {
0040         public:
0041         // Empty constructor,

0042         // defaults to compton scattering: gamma p -> gamma p

0043         reaction_kinematics()
0044         {
0045             _initial_state   = new two_body_state(0., M2_PROTON);
0046             _eps_gamma       = new polarization_vector(_initial_state);
0047             _target          = new dirac_spinor(_initial_state);
0048
0049             _final_state     = new two_body_state(0., M2_PROTON);
0050             _eps_vec         = new polarization_vector(_final_state);
0051             _recoil          = new dirac_spinor(_final_state);
0052         };
0053
0054         // Constructor with a set mX and JP

0055         // defaults to proton as baryon and real photon

0056         // string ID is deprecated but kept for backward compatibility

0057         reaction_kinematics(double mX, std::string id = "")
0058         : _mX(mX), _mX2(mX*mX)
0059         {
0060             _initial_state   = new two_body_state(0., M2_PROTON);
0061             _eps_gamma       = new polarization_vector(_initial_state);
0062             _target          = new dirac_spinor(_initial_state);
0063
0064             _final_state     = new two_body_state(mX*mX, M2_PROTON);
0065             _eps_vec         = new polarization_vector(_final_state);
0066             _recoil          = new dirac_spinor(_final_state);
0067         };
0068
0069
0070         // Constructor with a set mX and baryon mass mR

0071         // defaults to real photon

0072         reaction_kinematics(double mX, double mR)
0073         : _mX(mX), _mX2(mX*mX),
0074           _mR(mR), _mR2(mR*mR)
0075         {
0076             _initial_state   = new two_body_state(0., M2_PROTON);
0077             _eps_gamma       = new polarization_vector(_initial_state);
0078             _target          = new dirac_spinor(_initial_state);
0079
0080             _final_state     = new two_body_state(mX*mX, mR*mR);
0081             _eps_vec         = new polarization_vector(_final_state);
0082             _recoil          = new dirac_spinor(_final_state);
0083         };
0084
0085         // Constructor with a set mV and baryon mass mR

0086         // and massive incoming scatterer and target

0087         reaction_kinematics(double mX, double mR, double mT, double mB = 0.)
0088         : _mX(mX), _mX2(mX*mX),
0089           _mR(mR), _mR2(mR*mR),
0090           _mB(mB), _mB2(mB*mB),
0091           _mT(mT), _mT2(mT*mT)
0092         {
0093             _initial_state   = new two_body_state(mB*mB, mT*mT);
0094             _eps_gamma       = new polarization_vector(_initial_state);
0095             _target          = new dirac_spinor(_initial_state);
0096
0097             _final_state     = new two_body_state(mX*mX, mR*mR);
0098             _eps_vec         = new polarization_vector(_final_state);
0099             _recoil          = new dirac_spinor(_final_state);
0100         };
0101
0102         // destructor

0103         ~reaction_kinematics()
0104         {
0105             delete _initial_state;
0106             delete _final_state;
0107             delete _eps_gamma;
0108             delete _eps_vec;
0109             delete _target;
0110             delete _recoil;
0111         }
0112
0113         // ---------------------------------------------------------------------------

0114         // Masses

0115
0116         double _mB = 0., _mB2 = 0.;       // mass and mass squared of the "beam"

0117         double _mX = 0., _mX2 = 0.;       // mass and mass squared of the produced particle

0118
0119         double _mT = M_PROTON, _mT2 = M2_PROTON;  // mass of the target, assumed to be proton unless overriden

0120         double _mR = M_PROTON, _mR2 = M2_PROTON;  // mass of the recoil baryon, assumed to be proton unless overriden

0121
0122         inline double Wth(){ return (_mX + _mR); }; // square root of the threshold

0123         inline double sth(){ return Wth() * Wth(); }; // final state threshold

0124
0125         // Change the meson mass

0126         inline void set_mX(double m)
0127         {
0128             _mX  = m;
0129             _mX2 = m*m;
0130
0131             // also update the meson mass in two_body_state

0132             _final_state->set_mV2(m*m);
0133         };
0134
0135         inline void set_mX2(double m2)
0136         {
0137             _mX  = sqrt(m2);
0138             _mX2 = m2;
0139
0140             // also update the meson mass in two_body_state

0141             _final_state->set_mV2(m2);
0142         };
0143
0144         // Change virtuality of the photon

0145         // Q2 > 0

0146         inline void set_Q2(double q2)
0147         {
0148             if (q2 < 0) { std::cout << "Caution! set_Q2(x) requires x > 0! \n"; }
0149             _mB2 = -q2;
0150             _initial_state->set_mV2(-q2);
0151         };
0152
0153         // ---------------------------------------------------------------------------

0154         // Quantum numbers of produced meson.

0155         std::array<int,2> _jp{{1,1}};
0156         inline void set_JP(int J, int P)
0157         {
0158             _jp = {J, P};
0159             _helicities = get_helicities(J);
0160             _nAmps = _helicities.size();
0161         };
0162
0163         // Helicity configurations

0164         // Defaults to spin-1

0165         // Photon [0], Incoming Proton [1], Produced meson [2], Outgoing Proton [3]

0166         int _nAmps = 24;
0167         std::vector< std::array<int, 4> > _helicities = SPIN_ONE_HELICITIES;
0168
0169         //--------------------------------------------------------------------------

0170         two_body_state * _initial_state,  * _final_state;
0171         polarization_vector * _eps_vec, * _eps_gamma;
0172         dirac_spinor * _target, * _recoil;
0173
0174         // Get s-channel scattering angle from invariants

0175         inline double z_s(double s, double t)
0176         {
0177             std::complex<double> qdotqp = _initial_state->momentum(s) * _final_state->momentum(s);
0178             std::complex<double> E1E3   = _initial_state->energy_V(s) * _final_state->energy_V(s);
0179
0180             double result = t - _mX2 - _mB2 + 2.*abs(E1E3);
0181             result /= 2. * abs(qdotqp);
0182
0183             return result;
0184         };
0185
0186         // Scattering angle in the s-channel

0187         // Use TMath::ACos instead of std::acos because its safer at the end points

0188         inline double theta_s(double s, double t)
0189         {
0190             return TMath::ACos( z_s(s, t) );
0191         };
0192
0193         // Invariant variables

0194         inline double t_man(double s, double theta)
0195         {
0196             std::complex<double> qdotqp = _initial_state->momentum(s) * _final_state->momentum(s);
0197             std::complex<double> E1E3   = _initial_state->energy_V(s) * _final_state->energy_V(s);
0198
0199             return _mX2 + _mB2 - 2. * abs(E1E3) + 2. * abs(qdotqp) * cos(theta);
0200         };
0201
0202         inline double u_man(double s, double theta)
0203         {
0204             return _mX2 + _mB2 + _mT2 + _mR2 - s - t_man(s, theta);
0205         };
0206
0207         // Scattering angles in t and u channel frames

0208         inline std::complex<double> z_t(double s, double theta)
0209         {
0210             double t = t_man(s, theta);
0211             std::complex<double> p_t = sqrt(XR * Kallen(t, _mT2, _mR2)) / sqrt(XR * 4. * t);
0212             std::complex<double> q_t = sqrt(XR * Kallen(t, _mX2, _mB2)) / sqrt(XR * 4. * t);
0213
0214             std::complex<double> result;
0215             result = 2. * s + t - _mT2 - _mR2 - _mX2 - _mB2; // s - u

0216             result /= 4. * p_t * q_t;
0217
0218             return result;
0219         };
0220     };
0221 };
0222
0223 #endif