snippets/Exclusive/ePIC_ExcKinUtils.h

0001 //---------------------------------------------------------------------------------------
0002 //
0003 // Utility functions; calculation of kinematic quantities for exclusive ePIC analyses
0004 //
0005 // Author: O. Jevons, 27/02/25
0006 //
0007 //---------------------------------------------------------------------------------------
0008
0009 #include "TMath.h"
0010
0011 // Aliases for common 3/4-vector types
0012 using P3EVector=ROOT::Math::PxPyPzEVector;
0013 using P3MVector=ROOT::Math::PxPyPzMVector;
0014 using MomVector=ROOT::Math::DisplacementVector3D<ROOT::Math::Cartesian3D<Double_t>,ROOT::Math::DefaultCoordinateSystemTag>;
0015
0016 //-----------------------------------------------------------------------------------------------------------------------------
0017 // FUNCTION DEFINITIONS
0018 //
0019 // NOTE: Templating applied for brevity
0020 //       4-vector functions are valid for any type which contains the operators E(), P(), Pt() and M2()
0021 //       This includes TLorentzVector (legacy) and ALL variants of ROOT::Math::LorentzVector class
0022 //
0023 //-----------------------------------------------------------------------------------------------------------------------------
0024
0025
0026 //-----------------------------------------------------------------------------------------------------------------------------
0027 // FUNCTION DEFINITIONS
0028 //-----------------------------------------------------------------------------------------------------------------------------
0029
0030 // Calculate energy from momentum and mass
0031 // 1. Using vector structures for momentum
0032 // Works for ANY structure which contains Mag2() operator
0033 template<typename P>
0034 Double_t calcE(const P& mom, const Float_t& M){
0035   return TMath::Sqrt(mom.Mag2() + TMath::Power(M,2));
0036 }
0037 // 2. Using separate floats for momentum components
0038 Double_t calcE(const Float_t& px, const Float_t& py, const Float_t& pz, const Float_t& M){
0039   return TMath::Sqrt(TMath::Power(px,2) + TMath::Power(py,2) + TMath::Power(pz,2) + TMath::Power(M,2));
0040 }
0041
0042
0043 // Calculate Mandelstam t - BABE method using tRECO conventions
0044 // Uses incoming proton BEam and scattered BAryon 4-vectors
0045 // Another way of saying t = -(p' - p)^2
0046 //--------------------------------------------------------------
0047 // NEEDS: p, p'
0048 //--------------------------------------------------------------
0049 template <typename V>
0050 Double_t calcT_BABE(const V& be, const V& ba){
0051   double t = (ba - be).M2();
0052
0053   return TMath::Abs(t);
0054 }
0055
0056 // Calculate Mandelstam t - eX method using tRECO conventions
0057 // Uses difference between the beam and scattered electron and all of the (non-scattered) final state
0058 // e.g. for DVCS, X is a photon; for electroproduction, it is the species of interest; etc...
0059 //--------------------------------------------------------------
0060 // NEEDS: e, e', X   [q, X]
0061 //--------------------------------------------------------------
0062 // 1. Separate vectors for beam and scattered electrons
0063 template <typename V>
0064 Double_t calcT_eX(const V& e, const V& ep, const V& X){
0065   double t = (e - ep - X).M2();
0066   return TMath::Abs(t);
0067 }
0068 // 2. Giving virtual photon vector directly
0069 template <typename V>
0070 Double_t calcT_eX(const V& q, const V& X){
0071   double t = (q - X).M2();
0072   return TMath::Abs(t);
0073 }
0074
0075 // Calculate Mandelstam t - eXBA method using tRECO conventions
0076 // Include final state baryon information into eX method
0077 //--------------------------------------------------------------
0078 // NEEDS: e, e', p', X
0079 // CANNOT JUST GIVE q-VECTOR; NEED INFO. FROM e'
0080 //--------------------------------------------------------------
0081 template <typename V>
0082 Double_t calcT_eXBA(const V& e, const V& ep, const V& pp, const V& X){
0083   // Extract info. from vectors - e' energy and theta
0084   double E_ep = ep.E();
0085   double theta_ep = ep.Theta();
0086
0087   // Intermediate calculations - q vector and HFS sigma
0088   P3EVector q(e.X()-ep.X(), e.Y()-ep.Y(), e.Z()-ep.Z(), e.E()-ep.E());
0089   double sigma_h = (pp+X).E() - (pp+X).Z();
0090   double sigterm = sigma_h/2;
0091   double eterm = (E_ep*(1+TMath::Cos(theta_ep)))/2;
0092
0093   P3EVector pcorr(q.X(), q.Y(), -sigterm-eterm, sigterm-eterm);
0094
0095   double t = (pcorr - X).M2();
0096   return TMath::Abs(t);
0097 }
0098
0099 // Calculate Mandelstam t - eXBE method using tRECO conventions
0100 // Include initial beam proton information into eX method
0101 //--------------------------------------------------------------
0102 // NEEDS: e, p, ep, pp, X    [q, p, pp, X]
0103 //--------------------------------------------------------------
0104 // 1. Separate vectors for beam and scattered electrons
0105 template<typename V>
0106 Double_t calcT_eXBE(const V& e, const V& p, const V& ep, const V& pp, const V& X){
0107   // Calculate 'missing' momentum, ignoring scattered baryon vector
0108   P3EVector p4miss((e+p-ep-X).X(),(e+p-ep-X).Y(),(e+p-ep-X).Z(),(e+p-ep-X).E());
0109
0110   // Define corrected momentum vector using missing momentum and scattered baryon mass
0111   Float_t pmiss_mag = p4miss.Vect().R();
0112   Float_t pcorr_mag = TMath::Sqrt(TMath::Power(pmiss_mag,2) + TMath::Power(pp.M(),2));
0113   P3EVector pcorr(p4miss.Vect().X(), p4miss.Vect().Y(), p4miss.Vect().Z(), pcorr_mag);
0114
0115   double t = (pcorr-p).M2();
0116   return TMath::Abs(t);
0117 }
0118 // 2. Giving virtual photon vector directly
0119 template<typename V>
0120 Double_t calcT_eXBE(const V& p, const V& q, const V& pp, const V& X){
0121   // Calculate 'missing' momentum, ignoring scattered baryon vector
0122   P3EVector p4miss((p+q-X).X(),(p+q-X).Y(),(p+q-X).Z(),(p+q-X).E());
0123
0124   // Define corrected momentum vector using missing momentum and scattered baryon mass
0125   Float_t pmiss_mag = p4miss.Vect().R();
0126   Float_t pcorr_mag = TMath::Sqrt(TMath::Power(pmiss_mag,2) + TMath::Power(pp.M(),2));
0127   P3EVector pcorr(p4miss.Vect().X(), p4miss.Vect().Y(), p4miss.Vect().Z(), pcorr_mag);
0128
0129   double t = (pcorr-p).M2();
0130   return TMath::Abs(t);
0131 }
0132
0133 // Calculate Mandelstam t - eBABE method using tRECO conventions
0134 // Include electron (beam and scattered) information into BABE method
0135 //--------------------------------------------------------------
0136 // NEEDS: e, p, ep, pp    [q, p, pp]
0137 //--------------------------------------------------------------
0138 // 1. Separate vectors for beam and scattered electrons
0139 template<typename V>
0140 Double_t calcT_eBABE(const V& e, const V& p, const V& ep, const V& pp){
0141   // Calculate needed vectors
0142   P3EVector q(e.X()-ep.X(), e.Y()-ep.Y(), e.Z()-ep.Z(), e.E()-ep.E());
0143   P3EVector pcorr(-q.X(), -q.Y(), pp.Z(), pp.E());
0144
0145   double t = (pcorr - p).M2();
0146   return TMath::Abs(t);
0147 }
0148 // 2. Giving virtual photon vector directly
0149 template <typename V>
0150 Double_t calcT_eBABE(const V& p, const V& q, const V& pp){
0151   P3EVector pcorr(-q.X(), -q.Y(), pp.Z(), pp.E());
0152
0153   double t = (pcorr - p).M2();
0154   return TMath::Abs(t);
0155 }
0156
0157 // Calculate Mandelstam t - XBABE method using tRECO conventions
0158 // Includes further final state information into BABE method
0159 //--------------------------------------------------------------
0160 // NEEDS: p, pp, X
0161 // CANNOT PROVIDE HFS BY ITSELF
0162 //--------------------------------------------------------------
0163 template<typename V>
0164 Double_t calcT_XBABE(const V& p, const V& pp, const V& X){
0165   P3EVector pcorr(-X.X(), -X.Y(), pp.Z(), pp.E());
0166
0167   double t = (pcorr - p).M2();
0168   return TMath::Abs(t);
0169 }
0170
0171 // Calculate Mandelstam t - eXBABE method using tRECO conventions
0172 // Uses full event information
0173 //--------------------------------------------------------------
0174 // NEEDS: e, p, ep, pp, X    [q, p, pp, X]
0175 //--------------------------------------------------------------
0176 // 1. Separate vectors for beam and scattered electrons
0177 template<typename V>
0178 Double_t calcT_eXBABE(const V& e, const V& p, const V& ep, const V& pp, const V& X){
0179   // Calculate 'missing' momentum, ignoring scattered baryon vector
0180   P3EVector p4miss((e+p-ep-X).X(),(e+p-ep-X).Y(),(e+p-ep-X).Z(),(e+p-ep-X).E());
0181
0182   // Define corrected momentum vector using missing momentum and scattered baryon mass
0183   Float_t pmiss_mag = p4miss.Vect().R();
0184   ROOT::Math::Polar3DVector pcorr_vect(pmiss_mag, pp.Theta(), pp.Phi());
0185   Float_t pcorr_mag = TMath::Sqrt(TMath::Power(pmiss_mag,2) + TMath::Power(pp.M(),2));
0186
0187   P3EVector pcorr(pcorr_vect.X(), pcorr_vect.Y(), pcorr_vect.Z(), pcorr_mag);
0188
0189   double t = (pcorr-p).M2();
0190   return TMath::Abs(t);
0191 }
0192 // 2. Giving virtual photon vector directly
0193 template<typename V>
0194 Double_t calcT_eXBABE(const V& p, const V& q, const V& pp, const V& X){
0195   // Calculate 'missing' momentum, ignoring scattered baryon vector
0196   P3EVector p4miss((p+q-X).X(),(p+q-X).Y(),(p+q-X).Z(),(p+q-X).E());
0197
0198   // Define corrected momentum vector using missing momentum and scattered baryon mass
0199   Float_t pmiss_mag = p4miss.Vect().R();
0200   ROOT::Math::Polar3DVector pcorr_vect(pmiss_mag, pp.Theta(), pp.Phi());
0201   Float_t pcorr_mag = TMath::Sqrt(TMath::Power(pmiss_mag,2) + TMath::Power(pp.M(),2));
0202
0203   P3EVector pcorr(pcorr_vect.X(), pcorr_vect.Y(), pcorr_vect.Z(), pcorr_mag);
0204
0205   double t = (pcorr-p).M2();
0206   return TMath::Abs(t);
0207 }
0208
0209 // Calculate missing kinematics (mass/energy/momentum)
0210 // 3-body final state: ab->cdf
0211 // Missing momentum
0212 template <typename V>
0213 Double_t calcPMiss_3Body(const V& a, const V& b, const V& c, const V& d, const V& f){
0214   return (a+b-c-d-f).P();
0215 }
0216 // Missing transverse momentum
0217 template <typename V>
0218 Double_t calcPtMiss_3Body(const V& a, const V& b, const V& c, const V& d, const V& f){
0219   return (a+b-c-d-f).Pt();
0220 }
0221 // Missing energy
0222 template <typename V>
0223 Double_t calcEMiss_3Body(const V& a, const V& b, const V& c, const V& d, const V& f){
0224   return (a+b-c-d-f).E();
0225 }
0226 // Missing mass (squared)
0227 template <typename V>
0228 Double_t calcM2Miss_3Body(const V& a, const V& b, const V& c, const V& d, const V& f){
0229   Float_t fEMiss = (a+b-c-d-f).E();
0230   Float_t fPMiss = (a+b-c-d-f).P();
0231
0232   Float_t fM2Miss = TMath::Power(fEMiss,2) - TMath::Power(fPMiss,2);
0233   return fM2Miss;
0234 }
0235
0236 // 2-body final state: ab->cd
0237 // Missing momentum
0238 template <typename V>
0239 Double_t calcPMiss_2Body(const V& a, const V& b, const V& c, const V& d){
0240   return (a+b-c-d).P();
0241 }
0242 // Missing transverse momentum
0243 template <typename V>
0244 Double_t calcPtMiss_2Body(const V& a, const V& b, const V& c, const V& d){
0245   return (a+b-c-d).Pt();
0246 }
0247 // Missing energy
0248 template <typename V>
0249 Double_t calcEMiss_2Body(const V& a, const V& b, const V& c, const V& d){
0250   return (a+b-c-d).E();
0251 }
0252 // Missing mass (squared)
0253 template <typename V>
0254 Double_t calcM2Miss_2Body(const V& a, const V& b, const V& c, const V& d){
0255   Float_t fEMiss = (a+b-c-d).E();
0256   Float_t fPMiss = (a+b-c-d).P();
0257
0258   Float_t fM2Miss = TMath::Power(fEMiss,2) - TMath::Power(fPMiss,2);
0259   return fM2Miss;
0260 }