GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Born.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (c) 2003-2024, The GENIE Collaboration
4  For the full text of the license visit http://copyright.genie-mc.org
5 
6  Alfonso Garcia <aagarciasoto \at km3net.de>
7  IFIC & Harvard University
8 */
9 //____________________________________________________________________________
10 
11 #include "Physics/HELepton/XSection/Born.h"
12 #include "Framework/ParticleData/PDGLibrary.h"
13 #include "Framework/ParticleData/PDGCodes.h"
14 #include "Framework/Conventions/Constants.h"
15 
16 #include <TMath.h>
17 
18 #include <iostream>
19 
20 using namespace genie;
21 using namespace genie::constants;
22 
23 //____________________________________________________________________________
24 Born::Born()
25 {
26 
27  fGw = PDGLibrary::Instance()->Find(kPdgWM)->Width();
28  fGz = PDGLibrary::Instance()->Find(kPdgZ0)->Width();
29  fmw2c = TComplex(kMw2,fGw*kMw);
30  fmz2c = TComplex(kMz2,fGz*kMz);
31  TComplex rat = fmw2c/fmz2c;
32  fsw2 = TComplex(1.-rat.Re(),-rat.Im());
33  fcw2 = 1.-fsw2;
34  falpha = TMath::Sqrt(2.)*kGF/kPi * fmw2c * fsw2;
35 
36  fgae = -1./2. + 2.*fsw2;
37  fgbe = -1./2.;
38  fgav = 1./2.;
39 
40 }
41 //____________________________________________________________________________
42 Born::~Born()
43 {
44 
45 /*
46 Make sure the p3 is always the charged lepton:
47 
48 nu (p1) + lp (p2) -> lp (p3) + nu (p4)
49 
50  s-channel t-channel u-channel
51 1 \ / 3 1 --------- 3 1 ------\ / 3
52  ------ | | X
53 2 / \ 4 2 --------- 4 2 ------/ \ 4
54 */
55 
56 }
57 //____________________________________________________________________________
58 double Born::PXSecCCR(double s, double t, double mlin, double mlout)
59 {
60 
61  TComplex prop = falpha/fsw2/(s-fmw2c);
62 
63  return (t-mlout*mlout)*(t-mlin*mlin) * prop.Rho2();
64 
65 }
66 //____________________________________________________________________________
67 double Born::PXSecCCV(double s, double t, double mlin, double mlout)
68 {
69 
70  TComplex prop = falpha/fsw2/(t-fmw2c);
71 
72  return (s-mlout*mlout)*(s-mlin*mlin) * prop.Rho2();
73 
74 }
75 //____________________________________________________________________________
76 double Born::PXSecCCRNC(double s, double t, double mlin, double mlout)
77 {
78 
79  double u = GetU(mlin,mlout,s,t);
80 
81  TComplex a = fgav*(fgae-fgbe)/(u-fmz2c)/fcw2/fsw2;
82  TComplex b = fgav*(fgae+fgbe)/(u-fmz2c)/fcw2/fsw2 + 1./(s-fmw2c)/fsw2;
83  return falpha.Rho2() * ( (t-mlout*mlout)*(t-mlin*mlin)*b.Rho2() + (s-mlout*mlout)*(s-mlin*mlin)*a.Rho2() );
84 
85 }
86 //____________________________________________________________________________
87 double Born::PXSecCCVNC(double s, double t, double mlin, double mlout)
88 {
89 
90  double u = GetU(mlin,mlout,s,t);
91 
92  TComplex a = fgav*(fgae+fgbe)/(u-fmz2c)/fcw2/fsw2 + 1./(t-fmw2c)/fsw2;
93  TComplex b = fgav*(fgae-fgbe)/(u-fmz2c)/fcw2/fsw2;
94  return falpha.Rho2() * ( (t-mlout*mlout)*(t-mlin*mlin)*b.Rho2() + (s-mlout*mlout)*(s-mlin*mlin)*a.Rho2() );
95 
96 }
97 //____________________________________________________________________________
98 double Born::PXSecNCVnu(double s, double t, double mlin, double mlout)
99 {
100 
101  double u = GetU(mlin,mlout,s,t);
102 
103  TComplex a = fgav*(fgae+fgbe)/(u-fmz2c)/fcw2/fsw2;
104  TComplex b = fgav*(fgae-fgbe)/(u-fmz2c)/fcw2/fsw2;
105  return falpha.Rho2() * ( (t-mlout*mlout)*(t-mlin*mlin)*b.Rho2() + (s-mlout*mlout)*(s-mlin*mlin)*a.Rho2() );
106 
107 }
108 //____________________________________________________________________________
109 double Born::PXSecNCVnubar(double s, double t, double mlin, double mlout)
110 {
111 
112  double u = GetU(mlin,mlout,s,t);
113 
114  TComplex a = fgav*(fgae-fgbe)/(u-fmz2c)/fcw2/fsw2;
115  TComplex b = fgav*(fgae+fgbe)/(u-fmz2c)/fcw2;
116  return falpha.Rho2() * ( (t-mlout*mlout)*(t-mlin*mlin)*b.Rho2()/fsw2.Rho2() + (s-mlout*mlout)*(s-mlin*mlin)*a.Rho2() );
117 
118 }
119 //____________________________________________________________________________
120 double Born::PXSecPhoton(double s, double t, double mlout2)
121 {
122 
123  double u = kMw2 + mlout2 - s - t;
124 
125  double ME = 8*kPi2*falpha.Rho2()*s/kMw2/fsw2.Re()/TMath::Power(mlout2 - t,2)/TMath::Power(kMw2 - u,2)*
126  ( -2*(kMw2 - u)*TMath::Power(mlout2,3) - 2*TMath::Power(mlout2,2)*(-2*kMw2*u + u*(s + u) + TMath::Power(kMw2,2))
127  + mlout2*(-(kMw2*u*(4*s + 5*u)) + (s + u)*TMath::Power(kMw2,2) + 3*TMath::Power(kMw2,3) + (s + u)*TMath::Power(u,2))
128  + 2*kMw2*((3*s + u)*TMath::Power(kMw2,2) - TMath::Power(kMw2,3) + 4*u*TMath::Power(s,2) + 2*TMath::Power(s,3) + 3*s*TMath::Power(u,2)
129  + TMath::Power(u,3) - kMw2*TMath::Power(2*s + u,2)) );
130 
131  return TMath::Max(0.,ME);
132 
133 }
134 //____________________________________________________________________________
135 double Born::PXSecPhoton_T(double s12, double s13, double Q2, double ml2)
136 {
137  double ME2 = 0.0;
138  ME2 = (4*falpha.Rho2()*kPi2*(TMath::Power(ml2,4)*s12*(2*TMath::Power(kMw2,2)*TMath::Power(s12,2) - 2*kMw2*Q2*s12*s13 + TMath::Power(Q2,2)*s13*(-s12 + s13)) +
139  TMath::Power(ml2,3)*(-2*TMath::Power(kMw2,3)*TMath::Power(s12,3) + TMath::Power(Q2,2)*(s12 - s13)*s13*(-(Q2*s12) + TMath::Power(s12,2) + Q2*s13 - 3*s12*s13) +
140  2*TMath::Power(kMw2,2)*TMath::Power(s12,2)*(3*Q2*s12 - 2*TMath::Power(s12,2) + 2*Q2*s13 + s12*s13) +
141  kMw2*Q2*s12*(Q2*TMath::Power(s12,2) - Q2*TMath::Power(s13,2) - 2*s12*TMath::Power(s13,2))) +
142  TMath::Power(ml2,2)*(-6*TMath::Power(kMw2,4)*TMath::Power(s12,3) + 2*kMw2*Q2*
143  (-(Q2*s12*(s12 - 3*s13)) + TMath::Power(Q2,2)*(s12 - s13) + TMath::Power(s12,2)*(s12 - s13))*(s12 - s13)*s13 +
144  TMath::Power(Q2,2)*(s12 - s13)*TMath::Power(s13,2)*(-2*Q2*s12 + 2*TMath::Power(s12,2) + 2*Q2*s13 - 3*s12*s13) +
145  2*TMath::Power(kMw2,3)*TMath::Power(s12,2)*(2*TMath::Power(s12,2) - 3*Q2*(2*s12 + s13)) +
146  TMath::Power(kMw2,2)*s12*(2*TMath::Power(s12,2)*TMath::Power(s12 - s13,2) + TMath::Power(Q2,2)*(s12 - s13)*s13 -
147  2*Q2*s12*(TMath::Power(s12,2) - 8*s12*s13 - 2*TMath::Power(s13,2)))) -
148  2*TMath::Power(kMw2,2)*TMath::Power(s12,2)*(2*TMath::Power(kMw2,4)*s12 + TMath::Power(Q2,2)*TMath::Power(s13,2)*(-s12 + s13) -
149  2*TMath::Power(kMw2,3)*(2*TMath::Power(s12,2) - Q2*s13 + s12*s13) +
150  TMath::Power(kMw2,2)*(4*Q2*s12*(s12 - 2*s13) + TMath::Power(Q2,2)*(-s12 + s13) + 2*s12*TMath::Power(s12 + s13,2)) +
151  2*kMw2*s13*(TMath::Power(Q2,2)*(s12 - s13) - s12*(TMath::Power(s12,2) + TMath::Power(s13,2)) + Q2*(-TMath::Power(s12,2) + 2*s12*s13 + TMath::Power(s13,2)))) +
152  ml2*(10*TMath::Power(kMw2,5)*TMath::Power(s12,3) - TMath::Power(Q2,2)*(Q2 - s12)*TMath::Power(s12 - s13,2)*TMath::Power(s13,3) +
153  2*TMath::Power(kMw2,4)*TMath::Power(s12,2)*(3*Q2*s12 - 4*TMath::Power(s12,2) + 4*Q2*s13 - 3*s12*s13) +
154  kMw2*Q2*(s12 - s13)*TMath::Power(s13,2)*(2*TMath::Power(Q2,2)*(s12 - s13) + 2*TMath::Power(s12,2)*(s12 - s13) + Q2*s12*(-3*s12 + 5*s13)) +
155  TMath::Power(kMw2,3)*s12*(2*Q2*s12*(5*TMath::Power(s12,2) - 16*s12*s13 - TMath::Power(s13,2)) +
156  TMath::Power(Q2,2)*(-3*TMath::Power(s12,2) + 2*s12*s13 + TMath::Power(s13,2)) + 2*TMath::Power(s12,2)*(TMath::Power(s12,2) + 6*s12*s13 + TMath::Power(s13,2))) -
157  TMath::Power(kMw2,2)*s13*(TMath::Power(Q2,3)*TMath::Power(s12 - s13,2) - 2*TMath::Power(s12,3)*TMath::Power(s12 - s13,2) +
158  TMath::Power(Q2,2)*s12*(-5*TMath::Power(s12,2) + 8*s12*s13 - 3*TMath::Power(s13,2)) +
159  2*Q2*TMath::Power(s12,2)*(3*TMath::Power(s12,2) - 9*s12*s13 + 2*TMath::Power(s13,2))))))/(TMath::Power(kMw2,3)*TMath::Power(s12,2)*TMath::Power(ml2 - kMw2 + s13,2)*TMath::Power(ml2 - kMw2 - s12 + s13,2)*fsw2.Re());
160  return TMath::Max(0.,ME2);
161 }
162 //____________________________________________________________________________
163 double Born::PXSecPhoton_L(double s12, double s13, double Q2, double ml2)
164 {
165  double ME2 = 0.0;
166  ME2 = 2*falpha.Rho2()*Q2*TMath::Power(kMw2,-3)*kPi2*TMath::Power(s12,-2)*TMath::Power(ml2 - kMw2 + s13,-2)*TMath::Power(ml2 - kMw2 - s12 + s13,-2)*
167  ((s12 - s13)*TMath::Power(ml2,5)*TMath::Power(s12,2) - 2*s12*TMath::Power(ml2,4)*(2*kMw2*TMath::Power(s12,2) - (Q2 - s12)*(-3*s12*s13 + TMath::Power(s12,2) + 2*TMath::Power(s13,2))) +
168  TMath::Power(ml2,3)*(-2*(s12 - s13)*TMath::Power(kMw2,2)*TMath::Power(s12,2) +
169  2*kMw2*s12*(Q2*(9*s12*s13 - 5*TMath::Power(s12,2) - 2*TMath::Power(s13,2)) + s12*(-11*s12*s13 + 3*TMath::Power(s12,2) + 4*TMath::Power(s13,2))) +
170  (s12 - s13)*(TMath::Power(Q2,2)*TMath::Power(s12 - 2*s13,2) + TMath::Power(s12,2)*(-6*s12*s13 + TMath::Power(s12,2) + 6*TMath::Power(s13,2)) -
171  2*Q2*s12*(-5*s12*s13 + TMath::Power(s12,2) + 6*TMath::Power(s13,2)))) +
172  2*TMath::Power(ml2,2)*(4*(2*s12 - s13)*TMath::Power(kMw2,3)*TMath::Power(s12,2) +
173  (Q2 - s12)*s13*(Q2*(s12 - 2*s13) + s12*(-s12 + s13))*(-3*s12*s13 + TMath::Power(s12,2) + 2*TMath::Power(s13,2)) +
174  s12*TMath::Power(kMw2,2)*(Q2*(-11*s12*s13 + 7*TMath::Power(s12,2) - 2*TMath::Power(s13,2)) + s12*(15*s12*s13 - 11*TMath::Power(s12,2) + 2*TMath::Power(s13,2))) -
175  kMw2*((s12 - s13)*TMath::Power(Q2,2)*TMath::Power(s12 - 2*s13,2) + TMath::Power(s12,2)*(-11*s13*TMath::Power(s12,2) + TMath::Power(s12,3) + 20*s12*TMath::Power(s13,2) - 8*TMath::Power(s13,3)) -
176  2*Q2*s12*(-10*s13*TMath::Power(s12,2) + TMath::Power(s12,3) + 17*s12*TMath::Power(s13,2) - 6*TMath::Power(s13,3)))) +
177  4*TMath::Power(kMw2,2)*TMath::Power(s12,2)*(-(s13*(Q2 - 7*s12 + 4*s13)*TMath::Power(kMw2,2)) + (s12 - 2*s13)*TMath::Power(kMw2,3) + kMw2*(2*Q2 - s12 - 2*s13)*TMath::Power(s13,2) +
178  (-Q2 + s12)*TMath::Power(s13,3)) + ml2*(-15*(s12 - s13)*TMath::Power(kMw2,4)*TMath::Power(s12,2) +
179  2*s12*TMath::Power(kMw2,3)*(Q2*(7*s12*s13 - 3*TMath::Power(s12,2) + 2*TMath::Power(s13,2)) + s12*(-21*s12*s13 + 9*TMath::Power(s12,2) + 4*TMath::Power(s13,2))) +
180  TMath::Power(kMw2,2)*((s12 - s13)*TMath::Power(Q2,2)*TMath::Power(s12 - 2*s13,2) +
181  TMath::Power(s12,2)*(-31*s13*TMath::Power(s12,2) + TMath::Power(s12,3) + 60*s12*TMath::Power(s13,2) - 14*TMath::Power(s13,3)) -
182  2*Q2*s12*(-14*s13*TMath::Power(s12,2) + TMath::Power(s12,3) + 27*s12*TMath::Power(s13,2) - 6*TMath::Power(s13,3))) -
183  2*kMw2*s13*((s12 - s13)*TMath::Power(Q2,2)*TMath::Power(s12 - 2*s13,2) +
184  TMath::Power(s12,2)*(-8*s13*TMath::Power(s12,2) + TMath::Power(s12,3) + 11*s12*TMath::Power(s13,2) - 4*TMath::Power(s13,3)) +
185  Q2*s12*(15*s13*TMath::Power(s12,2) - 2*TMath::Power(s12,3) - 25*s12*TMath::Power(s13,2) + 10*TMath::Power(s13,3))) +
186  (s12 - s13)*TMath::Power(s13,2)*TMath::Power(Q2*(s12 - 2*s13) + s12*(-s12 + s13),2)))/fsw2.Re();
187  return TMath::Max(0.,ME2);
188 }
189 //____________________________________________________________________________
190 double Born::Lambda(double a, double b, double c)
191 {
192  return a*a + b*b + c*c - 2*a*b - 2*a*c - 2*b*c;
193 }
194 //____________________________________________________________________________
195 double Born::GetS(double mlin, double Enuin)
196 {
197  return 2. * mlin * Enuin + mlin*mlin;
198 }
199 //____________________________________________________________________________
200 double Born::GetT(double mlin, double mlout, double s, double costhCM)
201 {
202  //http://edu.itp.phys.ethz.ch/hs10/ppp1/PPP1_2.pdf [Sec 2.2.1]
203  double sum = mlin*mlin+mlout*mlout;
204  return ( (TMath::Sqrt(Lambda(s,0.,mlin*mlin)*Lambda(s,mlout*mlout,0.))*costhCM+mlin*mlin*mlout*mlout)/s + sum - s ) /2.;
205 }
206 //____________________________________________________________________________
207 double Born::GetU(double mlin, double mlout, double s, double t)
208 {
209  return mlin*mlin+mlout*mlout-s-t;
210 }
211 //____________________________________________________________________________
212 bool Born::IsInPhaseSpace(double mlin, double mlout, double Enuin, double Enuout)
213 {
214 
215  //https://arxiv.org/pdf/2007.14426.pdf [section 2.2]
216  double frac = Enuout/Enuin;
217  if ( frac < mlin/(mlin+2.*Enuin)+(mlout*mlout-mlin*mlin)/2./Enuin/(mlin+2.*Enuin) ) return false;
218  else if ( frac > 1.-(mlout*mlout-mlin*mlin)/2./Enuin/mlin ) return false;
219 
220  return true;
221 
222 }
223 
224 
225 
226 
genie::Born::GetT
double GetT(double mlin, double mlout, double s, double costhCM)
Definition: Born.cxx:200
genie::utils::kinematics::Q2
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:1077
genie::Born::PXSecPhoton_T
double PXSecPhoton_T(double s12, double s13, double Q2, double ml2)
Definition: Born.cxx:135
genie::kPdgWM
const int kPdgWM
Definition: PDGCodes.h:192
genie::units::s
static constexpr double s
Definition: Units.h:95
genie::Born::PXSecNCVnubar
double PXSecNCVnubar(double s, double t, double mlin, double mlout)
Definition: Born.cxx:109
genie::Born::IsInPhaseSpace
bool IsInPhaseSpace(double mlin, double mlout, double Enuin, double Enuout)
Definition: Born.cxx:212
genie::kPdgZ0
const int kPdgZ0
Definition: PDGCodes.h:190
genie::units::b
static constexpr double b
Definition: Units.h:78
genie::Born::PXSecCCVNC
double PXSecCCVNC(double s, double t, double mlin, double mlout)
Definition: Born.cxx:87
a
const double a
Definition: gUpMuFluxGen.cxx:164
genie::Born::Lambda
double Lambda(double a, double b, double c)
Definition: Born.cxx:190
genie::Born::PXSecPhoton_L
double PXSecPhoton_L(double s12, double s13, double Q2, double ml2)
Definition: Born.cxx:163
genie::Born::PXSecCCR
double PXSecCCR(double s, double t, double mlin, double mlout)
Definition: Born.cxx:58
genie::Born::GetS
double GetS(double mlin, double Enuin)
Definition: Born.cxx:195
genie::Born::PXSecCCRNC
double PXSecCCRNC(double s, double t, double mlin, double mlout)
Definition: Born.cxx:76
genie::PDGLibrary::Instance
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:68
genie::Born::PXSecCCV
double PXSecCCV(double s, double t, double mlin, double mlout)
Definition: Born.cxx:67
genie::Born::PXSecPhoton
double PXSecPhoton(double s, double t, double mlout2)
Definition: Born.cxx:120
genie::PDGLibrary::Find
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:86
genie::Born::PXSecNCVnu
double PXSecNCVnu(double s, double t, double mlin, double mlout)
Definition: Born.cxx:98
genie::Born::~Born
virtual ~Born()
Definition: Born.cxx:42
PDGCodes.h
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
genie::Born::GetU
double GetU(double mlin, double mlout, double s, double t)
Definition: Born.cxx:207
Generated on Wed Jul 17 2024 14:46:38 for GENIEGenerator by   doxygen 1.8.5