GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
PhysUtils.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  Costas Andreopoulos <c.andreopoulos \at cern.ch>
7  University of Liverpool
8 */
9 //____________________________________________________________________________
10 
11 #include <TMath.h>
12 #include <TVector3.h>
13 
16 
17 //___________________________________________________________________________
19  double m, const TLorentzVector & p4,
20  const TVector3 & p3hadr, double ct0 /*in fm*/, double K)
21 {
22 // m -> hadon mass (on-shell)
23 // p -> hadron momentum 4-vector (Lab)
24 // p3hadr -> hadronic-system momentum 3-vector (Lab)
25 
26  TVector3 p3 = p4.Vect(); // hadron's: p (px,py,pz)
27  double m2 = m*m; // m^2
28  double P = p4.P(); // |p|
29  double Pt = p3.Pt(p3hadr); // pT
30  double Pt2 = Pt*Pt; // pT^2
31  double fz = P*ct0*m/(m2+K*Pt2); // formation zone, in fm
32 
33  LOG("PhysUtil", pNOTICE)
34  << "Formation zone(|P| = " << P << " GeV, Pt = " << Pt
35  << " GeV, ct0 = " << ct0 << " fm, K = " << K << ") = " << fz << " fm";
36 
37  return fz;
38 }
39 //___________________________________________________________________________
40 double genie::utils::phys::R99118(double x, double Q2)
41 {
42 // Adapted from fortran code sent by V.Tvaskis
43 // PRL 98, 142301, 2007
44 //
45  double A[3] = { .06723, .46714, 1.89794 };
46  double B[3] = { .06347, .57468, -.35342 };
47  double C[3] = { .05992, .50885, 2.10807 };
48 
49  double consq2=2.; // Joining point of r1990 and re99118
50  double R=0;
51 
52  if(Q2 >= consq2) {
53  double FAC = 1+12.*(Q2/(1.+Q2))*(.125*.125/(x*x+.125*.125));
54  double RLOG = FAC/TMath::Log(Q2/.04);
55  double Q2thr = 5.*TMath::Power(1.-x,5);
56  double R_A = A[0]*RLOG + A[1] / TMath::Sqrt( TMath::Sqrt( TMath::Power(Q2,4)+TMath::Power(A[2],4) ));
57  double R_B = B[0]*RLOG + B[1]/Q2 + B[2]/( TMath::Power(Q2,2) + TMath::Power(.3,2) );
58  double R_C = C[0]*RLOG + C[1]/TMath::Sqrt( TMath::Power(Q2-Q2thr,2) + TMath::Power(C[2],2) );
59  R = (R_A+R_B+R_C)/3.;
60  }
61 
62  if(Q2 < consq2) {
63  double FAC = 1+12.*(consq2/(1.+consq2))*(.125*.125/(x*x+.125*.125));
64  double RLOG = FAC/TMath::Log(consq2/.04);
65  double Q2thr = 5.*TMath::Power(1.-x,5);
66  double R_A = A[0]*RLOG + A[1]/TMath::Sqrt(TMath::Sqrt( TMath::Power(consq2,4)+TMath::Power(A[2],4) ));
67  double R_B = B[0]*RLOG + B[1]/consq2 + B[2]/( TMath::Power(consq2,2) + TMath::Power(.3,2) );
68  double R_C = C[0]*RLOG + C[1]/TMath::Sqrt( TMath::Power(consq2-Q2thr,2) + TMath::Power(C[2],2) );
69  double norm=(R_A+R_B+R_C)/3.;
70  R=norm*(1.0-TMath::Exp(-Q2*9.212));
71  }
72  return R;
73 }
74 //___________________________________________________________________________
75 double genie::utils::phys::RWhitlow(double x, double Q2)
76 {
77 // Adapted from NeuGEN's rmodel_mod()
78 //
79 // Hugh's comments in original code:
80 // from NuTeV code provided by Donna Naples, May 2005
81 // added form for R below qsq=.35 GEV**2 from hep-ex/030807
82 //
83 // NuTEV comments:
84 //
85 //C Revised to make HTWIST select between more than two different
86 //C values of R.
87 //C
88 //C HTWIST = 'F' ==> RQCD - WITH LIMIT R > (2MX)**2/Q2
89 //C HTWIST = 'T' ==> RSLAC - WITH LIMIT R > (2MX)**2/Q2
90 //C HTWIST = '0' ==> R = 0 (NOT TRUE, BUT USEFUL FOR STUDIES)
91 //C HTWIST = '2' ==> R =.2 (NOT TRUE, BUT USEFUL FOR STUDIES)
92 //C HTWIST = 'C' ==> R-CALLAN-GROSS = (2MX)**2/Q2
93 //C
94 //C HTWIST = 'W' ==> R PARAMETERIZATION FROM WHITLOW'S THESIS
95 //C HTWIST = '+' ==> R PARAMETRIZATION FROM WHITLOW +15%
96 //C HTWIST = 'P' ==> R PARAMETRIZATION FROM WHITLOW +0.03
97 //C HTWIST = 'M' ==> R PARAMETRIZATION FROM WHITLOW -0.03
98 //C
99 //C 22-DEC-92 WGS
100 //
101  const double C2 = TMath::Power(0.125, 2);
102  const double B1 = 0.0635;
103  const double B2 = 0.5747;
104  const double B3 = -0.3534;
105 
106  double Q2R = TMath::Max(Q2, 0.35);
107  double ss = TMath::Log(Q2R/.04);
108 
109  double x2 = TMath::Power(x, 2.);
110  double Q4R = TMath::Power(Q2R, 2.);
111  double Q4 = TMath::Power(Q2, 2.);
112 
113  double theta = 1. + (12.*Q2R/(Q2R+1.)) * (C2/(C2+x2));
114  double R = (B1/ss)*theta + B2/Q2R + B3/(Q4R+.09);
115 
116  R = TMath::Max(R,0.);
117 
118  if(Q2 < 0.35) {
119  R *= ( 3.207*(Q2/(Q4+1.)) );
120  }
121  return R;
122 }
123 //___________________________________________________________________________
124 /*
125 void genie::utils::phys::ExtractStructFunc (
126  double x, double Q2, double d2sig_dxdy[3],
127  double & F1, double & F2, double & xF3)
128 {
129 // Solve the system of equations:
130 // (Sigma) = (A) x (SF) => SF = (A)^-1 x (Sigma)
131 //
132 
133  const double sign = 1;
134  const double M = kNucleonMass;
135 
136  TMatrixD A (3,3); // (row,col)
137  TMatrixD Sigma (3,1);
138 
139  A(0,0) = x * TMath::Power(y[0],2.);
140  A(1,0) = x * TMath::Power(y[1],2.);
141  A(2,0) = x * TMath::Power(y[2],2.);
142 
143  A(0,1) = 1. - y[0] - 0.5*M*x*y[0]/E[0];
144  A(1,1) = 1. - y[1] - 0.5*M*x*y[1]/E[1];
145  A(2,1) = 1. - y[2] - 0.5*M*x*y[2]/E[2];
146 
147  A(0,2) = sign * y[0] * (1.-0.5*y[0]);
148  A(1,2) = sign * y[1] * (1.-0.5*y[1]);
149  A(2,2) = sign * y[2] * (1.-0.5*y[2]);
150 
151  Sigma(0,0) = d2sig_dxdy[0] / (kGF2*M*E[0]/kPi);
152  Sigma(1,0) = d2sig_dxdy[1] / (kGF2*M*E[1]/kPi);
153  Sigma(2,0) = d2sig_dxdy[2] / (kGF2*M*E[2]/kPi);
154 
155  A.Invert();
156 
157  TMatrixD SF = A*Sigma;
158 
159  F1 = SF(0,0);
160  F2 = SF(1,0);
161  xF3 = SF(2,0);
162 }
163 //___________________________________________________________________________
164 */
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:1077
double FormationZone(double m, const TLorentzVector &p, const TVector3 &p3hadr, double ct0, double K)
Definition: PhysUtils.cxx:18
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
static constexpr double m2
Definition: Units.h:72
static constexpr double A
Definition: Units.h:74
double R99118(double x, double Q2)
PRL 98, 142301, 2007.
Definition: PhysUtils.cxx:40
#define pNOTICE
Definition: Messenger.h:61
double RWhitlow(double x, double Q2)
Definition: PhysUtils.cxx:75
static constexpr double m
Definition: Units.h:71