Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-10-25 09:06:15 

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 // Extra utilities - calculate 10ths of histogram contents
0031 template<typename h>
0032 void calcTenths(const h& hist){
0033   int iTenthBin{1};
0034 
0035   // Count over all bins in histogram
0036   for(int iHistBin{1}; iHistBin<hist->GetNbinsX(); iHistBin++){
0037     // Break out of loop after 90% mark (last 10% up to last bin)
0038     if(iTenthBin == 10) break;
0039 
0040     // Else calculate fraction of events up to ith bin
0041     float frac = (float)hist->Integral(1,iHistBin)/hist->Integral();
0042     
0043     // First time fraction of events is greater than the n*10% mark, print previous bin and increment n
0044     if(frac > (float)iTenthBin/10){
0045       cout<<(int)iTenthBin*10<<"% \t bin no. "<<iHistBin-1<<"\t bin centre "<<hist->GetBinCenter(iHistBin-1)<<endl;
0046       iTenthBin++;
0047     }
0048   }
0049 
0050   return;
0051 }
0052 
0053 
0054 // Calculate energy from momentum and mass
0055 // 1. Using vector structures for momentum
0056 // Works for ANY structure which contains Mag2() operator
0057 template<typename P>
0058 Double_t calcE(const P& mom, const Float_t& M){ 
0059   return TMath::Sqrt(mom.Mag2() + TMath::Power(M,2)); 
0060 }
0061 // 2. Using separate floats for momentum components
0062 Double_t calcE(const Float_t& px, const Float_t& py, const Float_t& pz, const Float_t& M){ 
0063   return TMath::Sqrt(TMath::Power(px,2) + TMath::Power(py,2) + TMath::Power(pz,2) + TMath::Power(M,2)); 
0064 }
0065 
0066 
0067 // Calculate Mandelstam t - BABE method using tRECO conventions
0068 // Uses incoming proton BEam and scattered BAryon 4-vectors
0069 // Another way of saying t = -(p' - p)^2
0070 //--------------------------------------------------------------
0071 // NEEDS: p, p'
0072 //--------------------------------------------------------------
0073 template <typename V>
0074 Double_t calcT_BABE(const V& be, const V& ba){
0075   double t = (ba - be).M2();
0076   
0077   return TMath::Abs(t);
0078 }
0079 
0080 // Calculate Mandelstam t - eX method using tRECO conventions
0081 // Uses difference between the beam and scattered electron and all of the (non-scattered) final state
0082 // e.g. for DVCS, X is a photon; for electroproduction, it is the species of interest; etc...
0083 //--------------------------------------------------------------
0084 // NEEDS: e, e', X   [q, X]
0085 //--------------------------------------------------------------
0086 // 1. Separate vectors for beam and scattered electrons
0087 template <typename V>
0088 Double_t calcT_eX(const V& e, const V& ep, const V& X){
0089   double t = (e - ep - X).M2();
0090   return TMath::Abs(t);
0091 }
0092 // 2. Giving virtual photon vector directly
0093 template <typename V>
0094 Double_t calcT_eX(const V& q, const V& X){
0095   double t = (q - X).M2();
0096   return TMath::Abs(t);
0097 }
0098 
0099 // Calculate Mandelstam t - eXBA method using tRECO conventions
0100 // Include final state baryon information into eX method
0101 //--------------------------------------------------------------
0102 // NEEDS: e, e', p', X 
0103 // CANNOT JUST GIVE q-VECTOR; NEED INFO. FROM e'
0104 //--------------------------------------------------------------
0105 template <typename V>
0106 Double_t calcT_eXBA(const V& e, const V& ep, const V& pp, const V& X){
0107   // Extract info. from vectors - e' energy and theta
0108   double E_ep = ep.E();
0109   double theta_ep = ep.Theta();
0110 
0111   // Intermediate calculations - q vector and HFS sigma
0112   P3EVector q(e.X()-ep.X(), e.Y()-ep.Y(), e.Z()-ep.Z(), e.E()-ep.E());
0113   double sigma_h = (pp+X).E() - (pp+X).Z();
0114   double sigterm = sigma_h/2;
0115   double eterm = (E_ep*(1+TMath::Cos(theta_ep)))/2;
0116   
0117   P3EVector pcorr(q.X(), q.Y(), -sigterm-eterm, sigterm-eterm);
0118   
0119   double t = (pcorr - X).M2();
0120   return TMath::Abs(t);
0121 }
0122 
0123 // Calculate Mandelstam t - eXBE method using tRECO conventions
0124 // Include initial beam proton information into eX method
0125 //--------------------------------------------------------------
0126 // NEEDS: e, p, ep, pp, X    [q, p, pp, X]
0127 //--------------------------------------------------------------
0128 // 1. Separate vectors for beam and scattered electrons
0129 template<typename V>
0130 Double_t calcT_eXBE(const V& e, const V& p, const V& ep, const V& pp, const V& X){
0131   // Calculate 'missing' momentum, ignoring scattered baryon vector
0132   P3EVector p4miss((e+p-ep-X).X(),(e+p-ep-X).Y(),(e+p-ep-X).Z(),(e+p-ep-X).E());
0133     
0134   // Define corrected momentum vector using missing momentum and scattered baryon mass
0135   Float_t pmiss_mag = p4miss.Vect().R();
0136   Float_t pcorr_mag = TMath::Sqrt(TMath::Power(pmiss_mag,2) + TMath::Power(pp.M(),2));
0137   P3EVector pcorr(p4miss.Vect().X(), p4miss.Vect().Y(), p4miss.Vect().Z(), pcorr_mag);
0138 
0139   double t = (pcorr-p).M2();
0140   return TMath::Abs(t);
0141 }
0142 // 2. Giving virtual photon vector directly
0143 template<typename V>
0144 Double_t calcT_eXBE(const V& p, const V& q, const V& pp, const V& X){
0145   // Calculate 'missing' momentum, ignoring scattered baryon vector
0146   P3EVector p4miss((p+q-X).X(),(p+q-X).Y(),(p+q-X).Z(),(p+q-X).E());
0147     
0148   // Define corrected momentum vector using missing momentum and scattered baryon mass
0149   Float_t pmiss_mag = p4miss.Vect().R();
0150   Float_t pcorr_mag = TMath::Sqrt(TMath::Power(pmiss_mag,2) + TMath::Power(pp.M(),2));
0151   P3EVector pcorr(p4miss.Vect().X(), p4miss.Vect().Y(), p4miss.Vect().Z(), pcorr_mag);
0152 
0153   double t = (pcorr-p).M2();
0154   return TMath::Abs(t);
0155 }
0156 // 3. Using scalar mass of scattered baryon (with e/e')
0157 template<typename V>
0158 Double_t calcT_eXBE(const V& e, const V& p, const V& ep, const Float_t& mb, const V& X){
0159   // Calculate 'missing' momentum, ignoring scattered baryon vector
0160   P3EVector p4miss((e+p-ep-X).X(),(e+p-ep-X).Y(),(e+p-ep-X).Z(),(e+p-ep-X).E());
0161     
0162   // Define corrected momentum vector using missing momentum and scattered baryon mass
0163   Float_t pmiss_mag = p4miss.Vect().R();
0164   Float_t pcorr_mag = TMath::Sqrt(TMath::Power(pmiss_mag,2) + TMath::Power(mb,2));
0165   P3EVector pcorr(p4miss.Vect().X(), p4miss.Vect().Y(), p4miss.Vect().Z(), pcorr_mag);
0166 
0167   double t = (pcorr-p).M2();
0168   return TMath::Abs(t);
0169 }
0170 // 4. Using scalar mass of scattered baryon (with q)
0171 template<typename V>
0172 Double_t calcT_eXBE(const V& p, const V& q, const Float_t& mb, const V& X){
0173   // Calculate 'missing' momentum, ignoring scattered baryon vector
0174   P3EVector p4miss((p+q-X).X(),(p+q-X).Y(),(p+q-X).Z(),(p+q-X).E());
0175     
0176   // Define corrected momentum vector using missing momentum and scattered baryon mass
0177   Float_t pmiss_mag = p4miss.Vect().R();
0178   Float_t pcorr_mag = TMath::Sqrt(TMath::Power(pmiss_mag,2) + TMath::Power(mb,2));
0179   P3EVector pcorr(p4miss.Vect().X(), p4miss.Vect().Y(), p4miss.Vect().Z(), pcorr_mag);
0180 
0181   double t = (pcorr-p).M2();
0182   return TMath::Abs(t);
0183 }
0184 
0185 // Calculate Mandelstam t - eBABE method using tRECO conventions
0186 // Include electron (beam and scattered) information into BABE method
0187 //--------------------------------------------------------------
0188 // NEEDS: e, p, ep, pp    [q, p, pp]
0189 //--------------------------------------------------------------
0190 // 1. Separate vectors for beam and scattered electrons
0191 template<typename V>
0192 Double_t calcT_eBABE(const V& e, const V& p, const V& ep, const V& pp){
0193   // Calculate needed vectors
0194   P3EVector q(e.X()-ep.X(), e.Y()-ep.Y(), e.Z()-ep.Z(), e.E()-ep.E());
0195   P3EVector pcorr(-q.X(), -q.Y(), pp.Z(), pp.E());
0196 
0197   double t = (pcorr - p).M2();
0198   return TMath::Abs(t);
0199 }
0200 // 2. Giving virtual photon vector directly
0201 template <typename V>
0202 Double_t calcT_eBABE(const V& p, const V& q, const V& pp){
0203   P3EVector pcorr(-q.X(), -q.Y(), pp.Z(), pp.E());
0204 
0205   double t = (pcorr - p).M2();
0206   return TMath::Abs(t);
0207 }
0208 
0209 // Calculate Mandelstam t - XBABE method using tRECO conventions
0210 // Includes further final state information into BABE method
0211 //--------------------------------------------------------------
0212 // NEEDS: p, pp, X
0213 // CANNOT PROVIDE HFS BY ITSELF
0214 //--------------------------------------------------------------
0215 template<typename V>
0216 Double_t calcT_XBABE(const V& p, const V& pp, const V& X){
0217   P3EVector pcorr(-X.X(), -X.Y(), pp.Z(), pp.E());
0218 
0219   double t = (pcorr - p).M2();
0220   return TMath::Abs(t);
0221 }
0222 
0223 // Calculate Mandelstam t - eXBABE method using tRECO conventions
0224 // Uses full event information
0225 //--------------------------------------------------------------
0226 // NEEDS: e, p, ep, pp, X    [q, p, pp, X]
0227 //--------------------------------------------------------------
0228 // 1. Separate vectors for beam and scattered electrons
0229 template<typename V>
0230 Double_t calcT_eXBABE(const V& e, const V& p, const V& ep, const V& pp, const V& X){
0231   // Calculate 'missing' momentum, ignoring scattered baryon vector
0232   P3EVector p4miss((e+p-ep-X).X(),(e+p-ep-X).Y(),(e+p-ep-X).Z(),(e+p-ep-X).E());
0233     
0234   // Define corrected momentum vector using missing momentum and scattered baryon mass
0235   Float_t pmiss_mag = p4miss.Vect().R();
0236   ROOT::Math::Polar3DVector pcorr_vect(pmiss_mag, pp.Theta(), pp.Phi());
0237   Float_t pcorr_mag = TMath::Sqrt(TMath::Power(pmiss_mag,2) + TMath::Power(pp.M(),2));
0238   
0239   P3EVector pcorr(pcorr_vect.X(), pcorr_vect.Y(), pcorr_vect.Z(), pcorr_mag);
0240 
0241   double t = (pcorr-p).M2();
0242   return TMath::Abs(t);
0243 }
0244 // 2. Giving virtual photon vector directly
0245 template<typename V>
0246 Double_t calcT_eXBABE(const V& p, const V& q, const V& pp, const V& X){
0247   // Calculate 'missing' momentum, ignoring scattered baryon vector
0248   P3EVector p4miss((p+q-X).X(),(p+q-X).Y(),(p+q-X).Z(),(p+q-X).E());
0249     
0250   // Define corrected momentum vector using missing momentum and scattered baryon mass
0251   Float_t pmiss_mag = p4miss.Vect().R();
0252   ROOT::Math::Polar3DVector pcorr_vect(pmiss_mag, pp.Theta(), pp.Phi());
0253   Float_t pcorr_mag = TMath::Sqrt(TMath::Power(pmiss_mag,2) + TMath::Power(pp.M(),2));
0254   
0255   P3EVector pcorr(pcorr_vect.X(), pcorr_vect.Y(), pcorr_vect.Z(), pcorr_mag);
0256 
0257   double t = (pcorr-p).M2();
0258   return TMath::Abs(t);
0259 }
0260 
0261 // Calculate Mandelstam t - eHe method from ECCE EDT paper
0262 // A. Bylinkin et al., Nucl. Instrum. Meth. A 1052, 168238 (2023); eq. 9
0263 // Needs full event information - give vectors in lab frame
0264 template<typename V>
0265 Double_t calcT_eHe(const V& e, const V& p, const V& ep, const V& pp, const V& X){
0266   // Extract intermediate quantities from 4-vectors
0267   double M = p.M();
0268   double nu = e.E() - ep.E();
0269   double Q2 = (e-ep).M2();
0270 
0271   // Calculate cos(theta between virtual and real photon)
0272   P3EVector q((e-ep).X(), (e-ep).Y(), (e-ep).Z(), (e-ep).E());
0273 
0274   double cTheta = TMath::Cos( (q-X).Theta() );
0275   
0276   double cosTerm = TMath::Sqrt((nu*nu)+Q2)*cTheta;
0277   double num = M*Q2 + (2*M*nu)*(nu-cosTerm);
0278   double den = M + nu - cosTerm;
0279 
0280   double t = num/den;
0281   return TMath::Abs(t);
0282 }
0283 
0284 // Calculate missing kinematics (mass/energy/momentum)
0285 // 3-body final state: ab->cdf
0286 // Missing momentum
0287 template <typename V>
0288 Double_t calcPMiss_3Body(const V& a, const V& b, const V& c, const V& d, const V& f){ 
0289   return (a+b-c-d-f).P(); 
0290 }
0291 // Missing transverse momentum
0292 template <typename V>
0293 Double_t calcPtMiss_3Body(const V& a, const V& b, const V& c, const V& d, const V& f){
0294   return (a+b-c-d-f).Pt(); 
0295 }
0296 // Missing energy
0297 template <typename V>
0298 Double_t calcEMiss_3Body(const V& a, const V& b, const V& c, const V& d, const V& f){
0299   return (a+b-c-d-f).E(); 
0300 }
0301 // Missing mass (squared)
0302 template <typename V>
0303 Double_t calcM2Miss_3Body(const V& a, const V& b, const V& c, const V& d, const V& f){
0304   Float_t fEMiss = (a+b-c-d-f).E();
0305   Float_t fPMiss = (a+b-c-d-f).P();
0306 
0307   Float_t fM2Miss = TMath::Power(fEMiss,2) - TMath::Power(fPMiss,2);
0308   return fM2Miss;
0309 }
0310 
0311 // 2-body final state: ab->cd
0312 // Missing momentum
0313 template <typename V>
0314 Double_t calcPMiss_2Body(const V& a, const V& b, const V& c, const V& d){
0315   return (a+b-c-d).P();
0316 }
0317 // Missing transverse momentum
0318 template <typename V>
0319 Double_t calcPtMiss_2Body(const V& a, const V& b, const V& c, const V& d){
0320   return (a+b-c-d).Pt();
0321 }
0322 // Missing energy
0323 template <typename V>
0324 Double_t calcEMiss_2Body(const V& a, const V& b, const V& c, const V& d){
0325   return (a+b-c-d).E();
0326 }
0327 // Missing mass (squared)
0328 template <typename V>
0329 Double_t calcM2Miss_2Body(const V& a, const V& b, const V& c, const V& d){
0330   Float_t fEMiss = (a+b-c-d).E();
0331   Float_t fPMiss = (a+b-c-d).P();
0332 
0333   Float_t fM2Miss = TMath::Power(fEMiss,2) - TMath::Power(fPMiss,2);
0334   return fM2Miss;
0335 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!