Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-04-02 09:14:12

0001 //---------------------------------------------------------------------------------------
0002 //
0003 // Utility functions; calculation of kinematic quantities for inclusive ePIC analyses
0004 //
0005 // Q2, x and y for all 5 methods used for the InclusiveKinematics branch of EICrecon
0006 //     - Electron, JB, DA, Sigma and eSigma
0007 //
0008 // Author: O. Jevons, 27/02/25
0009 //
0010 //---------------------------------------------------------------------------------------
0011 
0012 #include "TMath.h"
0013 
0014 // Aliases for common 3/4-vector types
0015 using P3EVector=ROOT::Math::PxPyPzEVector;
0016 using P3MVector=ROOT::Math::PxPyPzMVector;
0017 using MomVector=ROOT::Math::DisplacementVector3D<ROOT::Math::Cartesian3D<Double_t>,ROOT::Math::DefaultCoordinateSystemTag>;
0018 
0019 //-----------------------------------------------------------------------------------------------------------------------------
0020 // FUNCTION DEFINITIONS
0021 //
0022 // NOTE: Templating applied for brevity
0023 //       Functions are valid for any type which contains the operators X(), Y(), Z(), E(), Pt(), M2() and Dot()
0024 //       This includes TLorentzVector (legacy) and ALL variants of ROOT::Math::LorentzVector class
0025 //
0026 //
0027 // NOTE 2: Order of vectors (whichever are required) is ALWAYS as follows (for process ep -> e'p'X - change scattered baryon as necessary)
0028 //
0029 //         e -> p -> e' -> p' -> X
0030 //
0031 //-----------------------------------------------------------------------------------------------------------------------------
0032 
0033 // Calculate DIS kinematics - electron method
0034 // e + p -> e' + X
0035 // Q2
0036 template <typename V>
0037 Double_t calcQ2_Elec(const V& e, const V& ep){
0038   return -(e-ep).M2();
0039 }
0040 // x
0041 template <typename V>
0042 Double_t calcX_Elec(const V& e, const V& p, const V& ep){
0043   P3EVector q(e.X()-ep.X(), e.Y()-ep.Y(), e.Z()-ep.Z(), e.E()-ep.E());
0044   Double_t q2 = -q.Dot(q);
0045   Double_t den = 2*q.Dot(p);
0046 
0047   return q2/den;
0048 }
0049 // y
0050 template <typename V>
0051 Double_t calcY_Elec(const V& e, const V& p, const V& ep){
0052   P3EVector q(e.X()-ep.X(), e.Y()-ep.Y(), e.Z()-ep.Z(), e.E()-ep.E());
0053   Double_t num = p.Dot(q);
0054   Double_t den = p.Dot(e);
0055   
0056   return num/den;
0057 }
0058 // Calculate all 3 kinematic variables (Q2, x, y) simultaneously and pass them to holding variables
0059 template <typename V>
0060 void calcKin_Elec(const V& e, const V& p, const V& ep, Float_t &Q2, Float_t &x, Float_t &y){
0061   // Reset kinematic variables
0062   Q2 = 0;
0063   x = 0;
0064   y = 0;
0065 
0066   // Calculate kinematics
0067   P3EVector q(e.X()-ep.X(), e.Y()-ep.Y(), e.Z()-ep.Z(), e.E()-ep.E());
0068   Float_t Q2_e = -q.Dot(q);
0069   
0070   Float_t q_dot_p = q.Dot(p);
0071   Float_t x_e = Q2_e/(2*q_dot_p);
0072 
0073   Float_t e_dot_p = e.Dot(p);
0074   Float_t y_e = q_dot_p/e_dot_p;
0075 
0076   // Export variables
0077   Q2 = Q2_e;
0078   x = x_e;
0079   y = y_e;
0080 }
0081 
0082 // Calculate inclusive kinematics (Q2, x, y) with the JB method
0083 // SIDIS case ep -> e'p'X
0084 // Q2
0085 template <typename V>
0086 Double_t calcQ2_JB(const V& e, const V& pp, const V& X){
0087   // Intermediate variables
0088   Float_t pT2_h = TMath::Power((pp+X).X(),2) + TMath::Power((pp+X).Y(),2);
0089   Float_t sigma_h = (pp+X).E() - (pp+X).Z();
0090   
0091   // Kinematic variables
0092   Float_t y_jb = sigma_h / (2*e.E());
0093   Float_t Q2_jb = pT2_h / (1-y_jb);
0094 
0095   return Q2_jb;
0096 }
0097 // x
0098 template <typename V>
0099 Double_t calcX_JB(const V& e, const V& p, const V& pp, const V& X){
0100   // Intermediate variables
0101   Float_t pT2_h = TMath::Power((pp+X).X(),2) + TMath::Power((pp+X).Y(),2);
0102   Float_t sigma_h = (pp+X).E() - (pp+X).Z();
0103   
0104   // Kinematic variables
0105   Float_t y_jb = sigma_h / (2*e.E());
0106   Float_t Q2_jb = pT2_h / (1-y_jb);
0107   Float_t x_jb = Q2_jb / (4*e.E()*p.E()*y_jb);
0108   
0109   return x_jb;
0110 }
0111 // y
0112 template <typename V>
0113 Double_t calcY_JB(const V& e, const V& pp, const V& X){
0114   // Intermediate variables
0115   Float_t pT2_h = TMath::Power((pp+X).X(),2) + TMath::Power((pp+X).Y(),2);
0116   Float_t sigma_h = (pp+X).E() - (pp+X).Z();
0117 
0118   // Kinematic variables
0119   Float_t y_jb = sigma_h / (2*e.E());
0120 
0121   return y_jb;
0122 }
0123 // Calculate all 3 kinematic variables (Q2, x, y) simultaneously and pass them to holding variables 
0124 template <typename V>
0125 void calcKin_JB(const V& e, const V& p, const V& pp, const V& X, Float_t &Q2, Float_t &x, Float_t &y){
0126   // Reset kinematic variables
0127   Q2 = 0;
0128   x = 0;
0129   y = 0;
0130   
0131   // Intermediate variables
0132   Float_t pT2_h = TMath::Power((pp+X).X(),2) + TMath::Power((pp+X).Y(),2);
0133   Float_t sigma_h = (pp+X).E() - (pp+X).Z();
0134   
0135   // Kinematic variables
0136   Float_t y_jb = sigma_h / (2*e.E());
0137   Float_t Q2_jb = pT2_h / (1-y_jb);
0138   Float_t x_jb = Q2_jb / (4*e.E()*p.E()*y_jb);
0139   
0140   // Export kinematic variables
0141   Q2 = Q2_jb;
0142   x = x_jb;
0143   y = y_jb;
0144 }
0145 
0146 // DIS case ep -> e'X
0147 // Q2
0148 template <typename V>
0149 Double_t calcQ2_JB(const V& e, const V& HFS){
0150   // Intermediate variables
0151   Float_t pT2_h = TMath::Power(HFS.X(),2) + TMath::Power(HFS.Y(),2);
0152   Float_t sigma_h = HFS.E() - HFS.Z();
0153   
0154   // Kinematic variables
0155   Float_t y_jb = sigma_h / (2*e.E());
0156   Float_t Q2_jb = pT2_h / (1-y_jb);
0157 
0158   return Q2_jb;
0159 }
0160 // x
0161 template <typename V>
0162 Double_t calcX_JB(const V& e, const V& p, const V& HFS){
0163   // Intermediate variables
0164   Float_t pT2_h = TMath::Power(HFS.X(),2) + TMath::Power(HFS.Y(),2);
0165   Float_t sigma_h = HFS.E() - HFS.Z();
0166   
0167   // Kinematic variables
0168   Float_t y_jb = sigma_h / (2*e.E());
0169   Float_t Q2_jb = pT2_h / (1-y_jb);
0170   Float_t x_jb = Q2_jb / (4*e.E()*p.E()*y_jb);
0171   
0172   return x_jb;
0173 }
0174 // y
0175 template <typename V>
0176 Double_t calcY_JB(const V& e, const V& HFS){
0177   // Intermediate variables
0178   Float_t sigma_h = HFS.E() - HFS.Z();
0179   
0180   // Kinematic variables
0181   Float_t y_jb = sigma_h / (2*e.E());
0182 
0183   return y_jb;
0184 }
0185 // Calculate all 3 kinematic variables (Q2, x, y) simultaneously and pass them to holding variables 
0186 template <typename V>
0187 void calcKin_JB(const V& e, const V& p, const V& HFS, Float_t &Q2, Float_t &x, Float_t &y){
0188   // Reset kinematic variables
0189   Q2 = 0;
0190   x = 0;
0191   y = 0;
0192   
0193   // Intermediate variables
0194   Float_t pT2_h = TMath::Power(HFS.X(),2) + TMath::Power(HFS.Y(),2);
0195   Float_t sigma_h = HFS.E() - HFS.Z();
0196   
0197   // Kinematic variables
0198   Float_t y_jb = sigma_h / (2*e.E());
0199   Float_t Q2_jb = pT2_h / (1-y_jb);
0200   Float_t x_jb = Q2_jb / (4*e.E()*p.E()*y_jb);
0201   
0202   // Export kinematic variables
0203   Q2 = Q2_jb;
0204   x = x_jb;
0205   y = y_jb;
0206 }
0207 
0208 // Calculate inclusive kinematics (Q2, x, y) with the DA method
0209 // SIDIS case ep -> e'p'X
0210 // Q2
0211 template <typename V>
0212 Double_t calcQ2_DA(const V& e, const V& ep, const V& pp, const V& X){
0213   // Intermediate variables
0214   Float_t pT2_h = TMath::Power((pp+X).X(),2) + TMath::Power((pp+X).Y(),2);
0215   Float_t sigma_h = (pp+X).E() - (pp+X).Z();
0216   Float_t S2_h = sigma_h*sigma_h;
0217   Float_t gamma = TMath::ACos((pT2_h-S2_h)/(pT2_h+S2_h));
0218   
0219   // Kinematic variables
0220   Float_t Q2_da = 4*e.E()*e.E()*(1./TMath::Tan(ep.Theta()/2))*(1./(TMath::Tan(ep.Theta()/2) +  TMath::Tan(gamma/2)));
0221 
0222   return Q2_da;
0223 }
0224 // x
0225 template <typename V>
0226 Double_t calcX_DA(const V& e, const V& p, const V& ep, const V& pp, const V& X){
0227   // Intermediate variables
0228   Float_t pT2_h = TMath::Power((pp+X).X(),2) + TMath::Power((pp+X).Y(),2);
0229   Float_t sigma_h = (pp+X).E() - (pp+X).Z();
0230   Float_t S2_h = sigma_h*sigma_h;
0231   Float_t gamma = TMath::ACos((pT2_h-S2_h)/(pT2_h+S2_h));
0232   
0233   // Kinematic variables
0234   Float_t y_da = TMath::Tan(gamma/2) / (TMath::Tan(ep.Theta()/2) +  TMath::Tan(gamma/2));
0235   Float_t Q2_da = 4*e.E()*e.E()*(1./TMath::Tan(ep.Theta()/2))*(1./(TMath::Tan(ep.Theta()/2) +  TMath::Tan(gamma/2)));
0236   Float_t x_da = Q2_da/(4*e.E()*p.E()*y_da);
0237 
0238   return x_da;
0239 }
0240 // y
0241 template <typename V>
0242 Double_t calcY_DA(const V& ep, const V& pp, const V& X){
0243   // Intermediate variables
0244   Float_t pT2_h = TMath::Power((pp+X).X(),2) + TMath::Power((pp+X).Y(),2);
0245   Float_t sigma_h = (pp+X).E() - (pp+X).Z();
0246   Float_t S2_h = sigma_h*sigma_h;
0247   Float_t gamma = TMath::ACos((pT2_h-S2_h)/(pT2_h+S2_h));
0248   
0249   // Kinematic variables
0250   Float_t y_da = TMath::Tan(gamma/2) / (TMath::Tan(ep.Theta()/2) +  TMath::Tan(gamma/2));
0251 
0252   return y_da;
0253 }
0254 // Calculate all 3 kinematic variables (Q2, x, y) simultaneously and pass them to holding variables 
0255 template <typename V>
0256 void calcKin_DA(const V& e, const V& p, const V& ep, const V& pp, const V& X, Float_t &Q2, Float_t &x, Float_t &y){
0257   // Reset kinematic variables
0258   Q2 = 0;
0259   x = 0;
0260   y = 0;
0261  
0262   // Calculate intermediate quantities
0263   // Intermediate variables
0264   Float_t pT2_h = TMath::Power((pp+X).X(),2) + TMath::Power((pp+X).Y(),2);
0265   Float_t sigma_h = (pp+X).E() - (pp+X).Z();
0266   Float_t S2_h = sigma_h*sigma_h;
0267   Float_t gamma = TMath::ACos((pT2_h-S2_h)/(pT2_h+S2_h));
0268   
0269   // Kinematic variables
0270   Float_t y_da = TMath::Tan(gamma/2) / (TMath::Tan(ep.Theta()/2) +  TMath::Tan(gamma/2));
0271   Float_t Q2_da = 4*e.E()*e.E()*(1./TMath::Tan(ep.Theta()/2))*(1./(TMath::Tan(ep.Theta()/2) +  TMath::Tan(gamma/2)));
0272   Float_t x_da = Q2_da/(4*e.E()*p.E()*y_da);
0273 
0274   // Export kinematic variables
0275   Q2 = Q2_da;
0276   x = x_da;
0277   y = y_da;
0278 }
0279 
0280 // DIS case ep -> e'X
0281 // Q2
0282 template <typename V>
0283 Double_t calcQ2_DA(const V& e, const V& ep, const V& HFS){
0284   // Intermediate variables
0285   Float_t pT2_h = TMath::Power(HFS.X(),2) + TMath::Power(HFS.Y(),2);
0286   Float_t sigma_h = HFS.E() - HFS.Z();
0287   Float_t S2_h = sigma_h*sigma_h;
0288   Float_t gamma = TMath::ACos((pT2_h-S2_h)/(pT2_h+S2_h));
0289   
0290   // Kinematic variables
0291   Float_t Q2_da = 4*e.E()*e.E()*(1./TMath::Tan(ep.Theta()/2))*(1./(TMath::Tan(ep.Theta()/2) +  TMath::Tan(gamma/2)));
0292 
0293   return Q2_da;
0294 }
0295 // x
0296 template <typename V>
0297 Double_t calcX_DA(const V& e, const V& p, const V& ep, const V& HFS){
0298   // Intermediate variables
0299   Float_t pT2_h = TMath::Power(HFS.X(),2) + TMath::Power(HFS.Y(),2);
0300   Float_t sigma_h = HFS.E() - HFS.Z();
0301   Float_t S2_h = sigma_h*sigma_h;
0302   Float_t gamma = TMath::ACos((pT2_h-S2_h)/(pT2_h+S2_h));
0303   
0304   // Kinematic variables
0305   Float_t y_da = TMath::Tan(gamma/2) / (TMath::Tan(ep.Theta()/2) +  TMath::Tan(gamma/2));
0306   Float_t Q2_da = 4*e.E()*e.E()*(1./TMath::Tan(ep.Theta()/2))*(1./(TMath::Tan(ep.Theta()/2) +  TMath::Tan(gamma/2)));
0307   Float_t x_da = Q2_da/(4*e.E()*p.E()*y_da);
0308 
0309   return x_da;
0310 }
0311 // y
0312 template <typename V>
0313 Double_t calcY_DA(const V& ep, const V& HFS){
0314   // Intermediate variables
0315   Float_t pT2_h = TMath::Power(HFS.X(),2) + TMath::Power(HFS.Y(),2);
0316   Float_t sigma_h = HFS.E() - HFS.Z();
0317   Float_t S2_h = sigma_h*sigma_h;
0318   Float_t gamma = TMath::ACos((pT2_h-S2_h)/(pT2_h+S2_h));
0319   
0320   // Kinematic variables
0321   Float_t y_da = TMath::Tan(gamma/2) / (TMath::Tan(ep.Theta()/2) +  TMath::Tan(gamma/2));
0322 
0323   return y_da;
0324 }
0325 // Calculate all 3 kinematic variables (Q2, x, y) simultaneously and pass them to holding variables 
0326 template <typename V>
0327 void calcKin_DA(const V& e, const V& p, const V& ep, const V& HFS, Float_t &Q2, Float_t &x, Float_t &y){
0328   // Reset kinematic variables
0329   Q2 = 0;
0330   x = 0;
0331   y = 0;
0332  
0333   // Calculate intermediate quantities
0334   // Intermediate variables
0335   Float_t pT2_h = TMath::Power(HFS.X(),2) + TMath::Power(HFS.Y(),2);
0336   Float_t sigma_h = HFS.E() - HFS.Z();
0337   Float_t S2_h = sigma_h*sigma_h;
0338   Float_t gamma = TMath::ACos((pT2_h-S2_h)/(pT2_h+S2_h));
0339   
0340   // Kinematic variables
0341   Float_t y_da = TMath::Tan(gamma/2) / (TMath::Tan(ep.Theta()/2) +  TMath::Tan(gamma/2));
0342   Float_t Q2_da = 4*e.E()*e.E()*(1./TMath::Tan(ep.Theta()/2))*(1./(TMath::Tan(ep.Theta()/2) +  TMath::Tan(gamma/2)));
0343   Float_t x_da = Q2_da/(4*e.E()*p.E()*y_da);
0344 
0345   // Export kinematic variables
0346   Q2 = Q2_da;
0347   x = x_da;
0348   y = y_da;
0349 }
0350 
0351 // Calculate inclusive kinematics (Q2, x, y) with the Sigma method
0352 // SIDIS case ep -> e'p'X
0353 // Q2
0354 template <typename V>
0355 Double_t calcQ2_Sigma(const V& ep, const V& pp, const V& X){
0356   // Intermediate variables
0357   Float_t pT2_h = TMath::Power((pp+X).X(),2) + TMath::Power((pp+X).Y(),2);
0358   Float_t sigma_h = (pp+X).E() - (pp+X).Z();
0359   Float_t pT2_e = ep.Pt() * ep.Pt();
0360   Float_t sigma_e = ep.E() - ep.Z();
0361   Float_t sigma_tot = sigma_e + sigma_h;
0362 
0363   // Kinematic variables
0364   Float_t y_sigma = sigma_h/sigma_tot;
0365   Float_t Q2_sigma = pT2_e / (1-y_sigma);
0366   
0367   return Q2_sigma;
0368 }
0369 // x
0370 template <typename V>
0371 Double_t calcX_Sigma(const V& e, const V& p, const V& ep, const V& pp, const V& X){
0372   // Intermediate variables
0373   Float_t pT2_h = TMath::Power((pp+X).X(),2) + TMath::Power((pp+X).Y(),2);
0374   Float_t sigma_h = (pp+X).E() - (pp+X).Z();
0375   Float_t pT2_e = ep.Pt() * ep.Pt();
0376   Float_t sigma_e = ep.E() - ep.Z();
0377   Float_t sigma_tot = sigma_e + sigma_h;
0378 
0379   // Kinematic variables
0380   Float_t y_sigma = sigma_h/sigma_tot;
0381   Float_t Q2_sigma = pT2_e / (1-y_sigma);
0382   Float_t x_sigma = Q2_sigma / (4*e.E()*p.E()*y_sigma);
0383   
0384   return x_sigma;
0385 }
0386 // y
0387 template <typename V>
0388 Double_t calcY_Sigma(const V& ep, const V& pp, const V& X){
0389   // Intermediate variables
0390   Float_t pT2_h = TMath::Power((pp+X).X(),2) + TMath::Power((pp+X).Y(),2);
0391   Float_t sigma_h = (pp+X).E() - (pp+X).Z();
0392   Float_t sigma_e = ep.E() - ep.Z();
0393   Float_t sigma_tot = sigma_e + sigma_h;
0394 
0395   // Kinematic variables
0396   Float_t y_sigma = sigma_h/sigma_tot;
0397 
0398   return y_sigma;
0399 }
0400 // Calculate all 3 kinematic variables (Q2, x, y) simultaneously and pass them to holding variables 
0401 template <typename V>
0402 void calcKin_Sigma(const V& e, const V& p, const V& ep, const V& pp, const V& X, Float_t &Q2, Float_t &x, Float_t &y){
0403   // Reset kinematic variables
0404   Q2 = 0;
0405   x = 0;
0406   y = 0;
0407   
0408   // Intermediate variables
0409   Float_t pT2_h = TMath::Power((pp+X).X(),2) + TMath::Power((pp+X).Y(),2);
0410   Float_t sigma_h = (pp+X).E() - (pp+X).Z();
0411   Float_t pT2_e = ep.Pt() * ep.Pt();
0412   Float_t sigma_e = ep.E() - ep.Z();
0413   Float_t sigma_tot = sigma_e + sigma_h;
0414 
0415   // Kinematic variables
0416   Float_t y_sigma = sigma_h/sigma_tot;
0417   Float_t Q2_sigma = pT2_e / (1-y_sigma);
0418   Float_t x_sigma = Q2_sigma / (4*e.E()*p.E()*y_sigma);
0419   
0420   // Export kinematic variables
0421   Q2 = Q2_sigma;
0422   x = x_sigma;
0423   y = y_sigma;
0424 }
0425 
0426 // DIS case ep -> e'X
0427 // Q2
0428 template <typename V>
0429 Double_t calcQ2_Sigma(const V& ep, const V& HFS){
0430   // Intermediate variables
0431   Float_t pT2_h = TMath::Power(HFS.X(),2) + TMath::Power(HFS.Y(),2);
0432   Float_t sigma_h = HFS.E() - HFS.Z();
0433   Float_t pT2_e = ep.Pt() * ep.Pt();
0434   Float_t sigma_e = ep.E() - ep.Z();
0435   Float_t sigma_tot = sigma_e + sigma_h;
0436 
0437   // Kinematic variables
0438   Float_t y_sigma = sigma_h/sigma_tot;
0439   Float_t Q2_sigma = pT2_e / (1-y_sigma);
0440   
0441   return Q2_sigma;
0442 }
0443 // x
0444 template <typename V>
0445 Double_t calcX_Sigma(const V& e, const V& p, const V& ep, const V& HFS){
0446   // Intermediate variables
0447   Float_t pT2_h = TMath::Power(HFS.X(),2) + TMath::Power(HFS.Y(),2);
0448   Float_t sigma_h = HFS.E() - HFS.Z();
0449   Float_t pT2_e = ep.Pt() * ep.Pt();
0450   Float_t sigma_e = ep.E() - ep.Z();
0451   Float_t sigma_tot = sigma_e + sigma_h;
0452 
0453   // Kinematic variables
0454   Float_t y_sigma = sigma_h/sigma_tot;
0455   Float_t Q2_sigma = pT2_e / (1-y_sigma);
0456   Float_t x_sigma = Q2_sigma / (4*e.E()*p.E()*y_sigma);
0457   
0458   return x_sigma;
0459 }
0460 // y
0461 template <typename V>
0462 Double_t calcY_Sigma(const V& ep, const V& HFS){
0463   // Intermediate variables
0464   Float_t pT2_h = TMath::Power(HFS.X(),2) + TMath::Power(HFS.Y(),2);
0465   Float_t sigma_h = HFS.E() - HFS.Z();
0466   Float_t sigma_e = ep.E() - ep.Z();
0467   Float_t sigma_tot = sigma_e + sigma_h;
0468 
0469   // Kinematic variables
0470   Float_t y_sigma = sigma_h/sigma_tot;
0471 
0472   return y_sigma;
0473 }
0474 // Calculate all 3 kinematic variables (Q2, x, y) simultaneously and pass them to holding variables 
0475 template <typename V>
0476 void calcKin_Sigma(const V& e, const V& p, const V& ep, const V& HFS, Float_t &Q2, Float_t &x, Float_t &y){
0477   // Reset kinematic variables
0478   Q2 = 0;
0479   x = 0;
0480   y = 0;
0481   
0482   // Intermediate variables
0483   Float_t pT2_h = TMath::Power(HFS.X(),2) + TMath::Power(HFS.Y(),2);
0484   Float_t sigma_h = HFS.E() - HFS.Z();
0485   Float_t pT2_e = ep.Pt() * ep.Pt();
0486   Float_t sigma_e = ep.E() - ep.Z();
0487   Float_t sigma_tot = sigma_e + sigma_h;
0488 
0489   // Kinematic variables
0490   Float_t y_sigma = sigma_h/sigma_tot;
0491   Float_t Q2_sigma = pT2_e / (1-y_sigma);
0492   Float_t x_sigma = Q2_sigma / (4*e.E()*p.E()*y_sigma);
0493   
0494   // Export kinematic variables
0495   Q2 = Q2_sigma;
0496   x = x_sigma;
0497   y = y_sigma;
0498 }
0499 
0500 // Calculate inclusive kinematics (Q2, x, y) with the e-Sigma method
0501 // SIDIS case ep -> e'p'X
0502 // Q2
0503 template <typename V>
0504 Double_t calcQ2_ESigma(const V& e, const V& ep){
0505   // Intermediate variables
0506   Float_t pT2_e = ep.Pt() * ep.Pt();
0507   Float_t sigma_e = ep.E() - ep.Z();
0508 
0509   // Electron-based kinematics
0510   Float_t y_e = 1 - sigma_e / (2*e.E());
0511   Float_t Q2_e = pT2_e / (1-y_e);
0512 
0513   return Q2_e;
0514 }
0515 // x
0516 template <typename V>
0517 Double_t calcX_ESigma(const V& e, const V& p, const V& ep, const V& pp, const V& X){
0518   // Intermediate variables
0519   Float_t sigma_h = (pp+X).E() - (pp+X).Z();
0520   Float_t pT2_e = ep.Pt() * ep.Pt();
0521   Float_t sigma_e = ep.E() - ep.Z();
0522   Float_t sigma_tot = sigma_e + sigma_h;
0523 
0524   // Hadron-based kinematics
0525   Float_t y_sigma = sigma_h/sigma_tot;
0526   Float_t Q2_sigma = pT2_e / (1-y_sigma);
0527   Float_t x_sigma = Q2_sigma / (4*e.E()*p.E()*y_sigma);
0528 
0529   return x_sigma;
0530 }
0531 // y
0532 template <typename V>
0533 Double_t calcY_ESigma(const V& e, const V& p, const V& ep, const V& pp, const V& X){
0534   // Intermediate variables
0535   Float_t sigma_h = (pp+X).E() - (pp+X).Z();
0536   Float_t pT2_e = ep.Pt() * ep.Pt();
0537   Float_t sigma_e = ep.E() - ep.Z();
0538   Float_t sigma_tot = sigma_e + sigma_h;
0539 
0540   // Electron-based kinematics
0541   Float_t y_e = 1 - sigma_e / (2*e.E());
0542   Float_t Q2_e = pT2_e / (1-y_e);
0543 
0544   // Hadron-based kinematics
0545   Float_t y_sigma = sigma_h/sigma_tot;
0546   Float_t Q2_sigma = pT2_e / (1-y_sigma);
0547   Float_t x_sigma = Q2_sigma / (4*e.E()*p.E()*y_sigma);
0548 
0549   // Mixed kinematics
0550   Float_t y_esigma = Q2_e / (4*e.E()*p.E()*x_sigma);
0551 
0552   return y_esigma;
0553 }
0554 // Calculate all 3 kinematic variables (Q2, x, y) simultaneously and pass them to holding variables 
0555 template <typename V>
0556 void calcKin_ESigma(const V& e, const V& p, const V& ep, const V& pp, const V& X, Float_t &Q2, Float_t &x, Float_t &y){
0557   // Reset kinematic variables
0558   Q2 = 0;
0559   x = 0;
0560   y = 0;
0561   
0562   // Intermediate variables
0563   Float_t sigma_h = (pp+X).E() - (pp+X).Z();
0564   Float_t pT2_e = ep.Pt() * ep.Pt();
0565   Float_t sigma_e = ep.E() - ep.Z();
0566   Float_t sigma_tot = sigma_e + sigma_h;
0567 
0568   // Electron-based kinematics
0569   Float_t y_e = 1 - sigma_e / (2*e.E());
0570   Float_t Q2_e = pT2_e / (1-y_e);
0571 
0572   // Hadron-based kinematics
0573   Float_t y_sigma = sigma_h/sigma_tot;
0574   Float_t Q2_sigma = pT2_e / (1-y_sigma);
0575   Float_t x_sigma = Q2_sigma / (4*e.E()*p.E()*y_sigma);
0576 
0577   // Mixed kinematics
0578   Float_t y_esigma = Q2_e / (4*e.E()*p.E()*x_sigma);
0579 
0580   // Export kinematics
0581   Q2 = Q2_e;
0582   x = x_sigma;
0583   y = y_esigma;
0584 }
0585 
0586 // DIS case ep -> e'X
0587 // NO NEED TO REDEFINE Q2 - ONLY DEPENDS ON BEAM AND SCATTERED ELECTRON
0588 // x
0589 template <typename V>
0590 Double_t calcX_ESigma(const V& e, const V& p, const V& ep, const V& HFS){
0591   // Intermediate variables
0592   Float_t sigma_h = HFS.E() - HFS.Z();
0593   Float_t pT2_e = ep.Pt() * ep.Pt();
0594   Float_t sigma_e = ep.E() - ep.Z();
0595   Float_t sigma_tot = sigma_e + sigma_h;
0596 
0597   // Hadron-based kinematics
0598   Float_t y_sigma = sigma_h/sigma_tot;
0599   Float_t Q2_sigma = pT2_e / (1-y_sigma);
0600   Float_t x_sigma = Q2_sigma / (4*e.E()*p.E()*y_sigma);
0601 
0602   return x_sigma;
0603 }
0604 // y
0605 template <typename V>
0606 Double_t calcY_ESigma(const V& e, const V& p, const V& ep, const V& HFS){
0607   // Intermediate variables
0608   Float_t sigma_h = HFS.E() - HFS.Z();
0609   Float_t pT2_e = ep.Pt() * ep.Pt();
0610   Float_t sigma_e = ep.E() - ep.Z();
0611   Float_t sigma_tot = sigma_e + sigma_h;
0612 
0613   // Electron-based kinematics
0614   Float_t y_e = 1 - sigma_e / (2*e.E());
0615   Float_t Q2_e = pT2_e / (1-y_e);
0616 
0617   // Hadron-based kinematics
0618   Float_t y_sigma = sigma_h/sigma_tot;
0619   Float_t Q2_sigma = pT2_e / (1-y_sigma);
0620   Float_t x_sigma = Q2_sigma / (4*e.E()*p.E()*y_sigma);
0621 
0622   // Mixed kinematics
0623   Float_t y_esigma = Q2_e / (4*e.E()*p.E()*x_sigma);
0624 
0625   return y_esigma;
0626 }
0627 // Calculate all 3 kinematic variables (Q2, x, y) simultaneously and pass them to holding variables 
0628 template <typename V>
0629 void calcKin_ESigma(const V& e, const V& p, const V& ep, const V& HFS, Float_t &Q2, Float_t &x, Float_t &y){
0630    // Reset kinematic variables
0631   Q2 = 0;
0632   x = 0;
0633   y = 0;
0634   
0635   // Intermediate variables
0636   Float_t sigma_h = HFS.E() - HFS.Z();
0637   Float_t pT2_e = ep.Pt() * ep.Pt();
0638   Float_t sigma_e = ep.E() - ep.Z();
0639   Float_t sigma_tot = sigma_e + sigma_h;
0640 
0641   // Electron-based kinematics
0642   Float_t y_e = 1 - sigma_e / (2*e.E());
0643   Float_t Q2_e = pT2_e / (1-y_e);
0644 
0645   // Hadron-based kinematics
0646   Float_t y_sigma = sigma_h/sigma_tot;
0647   Float_t Q2_sigma = pT2_e / (1-y_sigma);
0648   Float_t x_sigma = Q2_sigma / (4*e.E()*p.E()*y_sigma);
0649 
0650   // Mixed kinematics
0651   Float_t y_esigma = Q2_e / (4*e.E()*p.E()*x_sigma);
0652 
0653   // Export kinematics
0654   Q2 = Q2_e;
0655   x = x_sigma;
0656   y = y_esigma;
0657 }