GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
HadXSUtils.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 
16 
17 using namespace genie::constants;
18 
19 //____________________________________________________________________________
21  bool isChargedPion)
22 {
23 // Returns the interpolated inelastic pion-nucleon cross section.
24 // C++ adaptation of Hugh Gallagher's NeuGEN inel() function
25 
26  double mpi = kPionMass;
27  double mpi2 = kPionMass2;
28  if (!isChargedPion) {
29  mpi = kPi0Mass;
30  mpi2 = mpi * mpi;
31  }
32  double Epion2 = TMath::Power(Epion,2);
33  double P = TMath::Sqrt( TMath::Max(0.,Epion2-mpi2) );
34 
35  if(P<=0) return 0;
36 
37  double log10P = TMath::Log10(P);
38  int N = (int) ((log10P - kInelMinLog10P)/kIneldLog10P) + 1;
39 
40  double xs=0.;
41  if ((log10P - kInelMinLog10P) < 0.0) xs = (P/0.1059)*kInelSig[0];
42  else if (N>kInelNDataPoints-2) xs = kInelSig[kInelNDataPoints-1];
43  else {
44  double log10Pn = kInelMinLog10P + (N-1) * kIneldLog10P;
45  double delta = (kInelSig[N]-kInelSig[N-1])/kIneldLog10P;
46  xs = kInelSig[N-1] + delta * (log10P-log10Pn);
47  }
48  return (xs * units::mb);
49 }
50 //____________________________________________________________________________
52  bool isChargedPion)
53 {
54 // Returns the interpolated total pion-nucleon cross section.
55 // C++ adaptation of Hugh Gallagher's NeuGEN total() function.
56 
57  double mpi = kPionMass;
58  double mpi2 = kPionMass2;
59  if (!isChargedPion) {
60  mpi = kPi0Mass;
61  mpi2 = mpi * mpi;
62  }
63  double Epion2 = TMath::Power(Epion,2);
64  double P = TMath::Sqrt( TMath::Max(0.,Epion2-mpi2) );
65 
66  if(P<=0) return 0;
67 
68  double log10P = TMath::Log10(P);
69  int N = (int) ((log10P - kInelMinLog10P)/kIneldLog10P) + 1;
70 
71  double xs=0.;
72  if ((log10P - kInelMinLog10P) < 0.0) xs = (P/0.1059)*kTotSig[0];
73  else if (N>kInelNDataPoints-2) xs = kTotSig[kInelNDataPoints-1];
74  else {
75  double log10Pn = kTotMinLog10P + (N-1) * kTotdLog10P;
76  double delta = (kTotSig[N]-kTotSig[N-1])/kTotdLog10P;
77  xs = kTotSig[N-1] + delta * (log10P-log10Pn);
78  }
79  return (xs * units::mb);
80 }
81 //____________________________________________________________________________
83  bool isChargedPion)
84 {
85  const double total = PionNucleonXSec(Epion, true, isChargedPion);
86  const double elastic = PionNucleonXSec(Epion, false, isChargedPion);
87  return (total - elastic);
88 }
89 //____________________________________________________________________________
91  bool isChargedPion)
92 {
93  return PionNucleonXSec(Epion, true, isChargedPion);
94 }
95 //____________________________________________________________________________
96 double genie::utils::hadxs::berger::PionNucleonXSec(double Epion, bool get_total,
97  bool isChargedPion)
98 {
99  // Convert inputs from Genie to those expected by Berger's code:
100 
101  double mpi = kPionMass;
102  double mpi2 = kPionMass2;
103  if (!isChargedPion) {
104  mpi = kPi0Mass;
105  mpi2 = mpi * mpi;
106  }
107 
108  double Epion2 = TMath::Power(Epion,2);
109  double ppi = TMath::Sqrt( TMath::Max(0., Epion2 - mpi2) );
110 
111  if( ppi <= 0.0 ) return 0.0;
112 
113  int out;
114  if( get_total ) out = 0; // Total pion-nucleon cross-section
115  else out = 1; // Elastic pion-nucleon cross-section
116 
117  const double M_pi = mpi; // TODO: used to be kPionMass prior to charge checks
118  const double M_p = kProtonMass; // TODO: should be kNucleonMass ??
119  const double pi = kPi;
120 
121  // Now this is the Berger's code...
122 
123  double afit=1.0;
124  double afit2=1.9;
125  double afit3=0.27;
126  double afit4=0.34;
127  double afit5=0.75;
128  double afit6=1.7;
129 
130  double epi = TMath::Sqrt(M_pi*M_pi + ppi*ppi);
131  double sx = M_p*M_p + M_pi*M_pi + 2.0 * epi * M_p;
132  double Wx = TMath::Sqrt(sx);
133 
134  double s12x = TMath::Sqrt((sx - TMath::Power((M_pi+M_p),2)) * (sx - (M_pi-M_p)*(M_pi-M_p)));
135  double ppistx = s12x/2.0/Wx;
136 
137  double Wdel = 1.232;
138  double arg1 = (Wdel*Wdel - (M_p + M_pi)*(M_p + M_pi))*(Wdel*Wdel - (M_p -M_pi)*(M_p -M_pi));
139  if(arg1<=0) arg1 = 0.0;
140  double gamx = 0.15 * TMath::Power((6.0*ppistx),3) / (1.0 + (6*ppistx)*(6*ppistx));
141  double bwnrx = 1.0 / (4.0*(Wx - Wdel)*(Wx - Wdel) + gamx*gamx);
142  double f1x = afit * 0.3892 * 8.0 * pi / (ppistx*ppistx) * gamx*gamx * bwnrx;
143  double f4x = afit4 * 0.3892 * 8.0 * pi / (ppistx*ppistx) * gamx*gamx * bwnrx;
144 
145  double Wres2 = 1.98;
146  double Gres2 = 0.6;
147  double arg2 = (Wres2*Wres2 - (M_p+M_pi)*(M_p+M_pi)) * (Wres2*Wres2 - (M_p - M_pi)*(M_p - M_pi));
148  if(arg2<=0) arg2 = 0.0;
149  double pp2 = TMath::Sqrt(arg2)/2.0/Wres2;
150  double gam2x = Gres2 * TMath::Power((ppistx/pp2),3);
151  double bwnr2x = 1.0 / (4.0 * (Wx-Wres2)*(Wx-Wres2) + gam2x*gam2x);
152  double f2x = afit2 * 0.3892 * 8.0 * pi / (ppistx*ppistx) * gam2x*gam2x * bwnr2x;
153 
154  double Wres3 = 1.7;
155  double Gres3 = 0.27;
156  double arg3 = (Wres3*Wres3 - (M_p+M_pi)*(M_p+M_pi)) * (Wres3*Wres3 - (M_p - M_pi)*(M_p - M_pi));
157  if(arg3<=0) arg3 = 0.0;
158  double pp3 = TMath::Sqrt(arg3)/2.0/Wres3;
159  double gam3x = Gres3 * TMath::Power((ppistx/pp3),3);
160  double bwnr3x = 1.0 / (4.0 * (Wx-Wres3)*(Wx-Wres3) + gam3x*gam3x);
161  double f3x = afit3 * 0.3892 * 8.0 * pi / (ppistx*ppistx) * gam3x*gam3x * bwnr3x;
162 
163  double Wres5 = 1.52;
164  double Gres5 = 0.10;
165  double arg5 = (Wres5*Wres5 - (M_p+M_pi)*(M_p+M_pi)) * (Wres5*Wres5 - (M_p - M_pi)*(M_p - M_pi));
166  if(arg5<=0) arg5 = 0.0;
167  double pp5 = TMath::Sqrt(arg5)/2.0/Wres5;
168  double gam5x = Gres5 * (ppistx/pp5);
169  double bwnr5x = 1.0 / (4.0 * (Wx-Wres5)*(Wx-Wres5) + gam5x*gam5x);
170  double f5x = afit5 * 0.3892 * 8.0 * pi / (ppistx*ppistx) * gam5x*gam5x * bwnr5x;
171 
172  double Wres6 = 1.69;
173  double Gres6 = 0.140;
174  double arg6 = (Wres6*Wres6 - (M_p+M_pi)*(M_p+M_pi)) * (Wres6*Wres6 - (M_p - M_pi)*(M_p - M_pi));
175  if(arg6<=0) arg6 = 0.0;
176  double pp6 = TMath::Sqrt(arg6)/2.0/Wres6;
177  double gam6x = Gres6 * (ppistx/pp6);
178  double bwnr6x = 1.0 / (4.0 * (Wx-Wres6)*(Wx-Wres6) + gam6x*gam6x);
179  double f6x = afit6 * 0.3892 * 8.0 * pi / (ppistx*ppistx) * gam6x*gam6x * bwnr6x;
180 
181  double output = -9999.99;
182  if(out==0){
183  double bgplusx = 0.834 * ppi + 1.77 * ppi*ppi;
184  double sigplusx = f1x + 0.86 * f2x + 0.7 * f3x + bgplusx;
185  if(ppi>1.9) sigplusx = 20.1 + 15.4 / TMath::Sqrt(ppi);
186  double bgminx = 17.8 * ppi + 6.22 * ppi*ppi - 3.1 * ppi*ppi*ppi;
187  double sigminx = 0.94 * f4x + 0.65 *f5x + 0.61 * f6x + bgminx;
188  if(ppi>1.2) sigminx = 21.64 + 17.29 / TMath::Sqrt(ppi);
189  double sigma_t = (sigplusx + sigminx) / 2.0;
190  output = sigma_t;
191  }else if(out==1){
192  double sgpelx = f1x + 0.411 * f2x;
193  if(ppi>1.9) sgpelx = 2.38 + 17.55 / ppi;
194  double sgmelx = 0.275 * f4x + 0.268 * f5x + 0.281 * f6x + 11.83 * ppi - 3.39 * ppi*ppi;
195  if(ppi>1.5) sgmelx = 3.4 + 11.35 / ppi;
196  double sigma_el = (sgpelx + sgmelx) / 2.0;
197  output = sigma_el;
198  }
199 
200  // Berger's code over, convert to Genie units and return
201  return (output * units::mb);
202 }
203 //____________________________________________________________________________
204 // Berger code for calculating the pion-Carbon xsec. If A is not 12 the restults are extrapolated to a nucleus of size A.
205 int genie::utils::hadxs::berger::PionNucleusXSec(double tpi, double ppistar, double t_new, double A, double &tpilow, double &siglow, double &tpihigh, double &sighigh){
206 
207  //Berger code for Tpi<=1.0
208  //Returns the entire dsigma(pi + N -> pi + N)/dt term based on pi-Carbon scattering data
209  double binedges[13] = {0.000, 0.076, 0.080, 0.100, 0.148, 0.162, 0.226, 0.486, 0.584, 0.662, 0.766, 0.870, 1.000};
210  double parones[13] = {0.0, 11600.0, 14700.0, 18300.0, 21300.0, 22400.0, 16400.0, 5730.0, 4610.0, 4570.0, 4930.0, 5140.0, 5140.0};
211  double partwos[13] = {0.0, 116.0, 109.0, 89.8, 91.0, 89.2, 80.8, 54.6, 55.2, 58.4, 60.5, 62.2, 62.2};
212  double factors[13] = {0.0, 0.1612, 0.1662, 0.1906, 0.2452, 0.2604, 0.3273, 0.5784, 0.6682, 0.7384, 0.8304, 0.9206, 0.9206};
213 
214  if(tpi>binedges[12]) return 1;
215  int btu = 1;
216  while(true){
217  if(tpi>binedges[btu]) btu++;
218  else break;
219  }
220  btu--;
221  if(btu<0) btu = 0;
222 
223  tpilow = binedges[btu];
224  tpihigh = binedges[btu + 1];
225 
226  double dsigdzlow = -9999.99;
227  if(btu==0) dsigdzlow = 0.0;
228  else dsigdzlow = 2.0 * factors[btu]*factors[btu]
229  * parones[btu]*TMath::Power(A/12.0,1.3333333)
230  * TMath::Exp( -1.0 * partwos[btu]*TMath::Power(A/12.0,0.6666666) * t_new * factors[btu]*factors[btu] / (ppistar*ppistar) );
231  siglow = dsigdzlow;
232 
233  double dsigdzhigh = -9999.99;
234  btu++;
235  if(btu>12) btu = 12;
236  if(btu==0) dsigdzhigh = 0.0;
237  else dsigdzhigh = 2.0 * factors[btu]*factors[btu]
238  * parones[btu]*TMath::Power(A/12.0,1.3333333)
239  * TMath::Exp( -1.0 * partwos[btu]*TMath::Power(A/12.0,0.6666666) * t_new * factors[btu]*factors[btu] / (ppistar*ppistar) );
240  sighigh = dsigdzhigh;
241 
242  return 0;
243 }
244 //_____________________________________________________________________________
static const double kTotdLog10P
Definition: HadXSUtils.h:53
double TotalPionNucleonXSec(double Epion, bool isChargedPion=true)
Definition: HadXSUtils.cxx:90
double PionNucleonXSec(double Epion, bool get_total, bool isChargedPion=true)
Definition: HadXSUtils.cxx:96
static const double kTotMinLog10P
Definition: HadXSUtils.h:52
static const double kInelMinLog10P
Definition: HadXSUtils.h:34
static constexpr double A
Definition: Units.h:74
static const double kIneldLog10P
Definition: HadXSUtils.h:35
static constexpr double mb
Definition: Units.h:79
static const int kInelNDataPoints
Definition: HadXSUtils.h:33
double InelasticPionNucleonXSec(double Epion, bool isChargedPion=true)
Definition: HadXSUtils.cxx:20
static const double kTotSig[kTotNDataPoints]
Definition: HadXSUtils.h:54
double InelasticPionNucleonXSec(double Epion, bool isChargedPion=true)
Definition: HadXSUtils.cxx:82
static const double kInelSig[kInelNDataPoints]
Definition: HadXSUtils.h:36
double TotalPionNucleonXSec(double Epion, bool isChargedPion=true)
Definition: HadXSUtils.cxx:51
int PionNucleusXSec(double tpi, double ppistar, double t_new, double A, double &tpilow, double &siglow, double &tpihigh, double &sighigh)
Definition: HadXSUtils.cxx:205