Functions
genie::utils::hadxs::berger Namespace Reference

Functions

double InelasticPionNucleonXSec (double Epion, bool isChargedPion=true)
 
double TotalPionNucleonXSec (double Epion, bool isChargedPion=true)
 
double PionNucleonXSec (double Epion, bool get_total, bool isChargedPion=true)
 
int PionNucleusXSec (double tpi, double ppistar, double t_new, double A, double &tpilow, double &siglow, double &tpihigh, double &sighigh)
 

Function Documentation

double genie::utils::hadxs::berger::InelasticPionNucleonXSec ( double  Epion,
bool  isChargedPion = true 
)

Definition at line 82 of file HadXSUtils.cxx.

84 {
85  const double total = PionNucleonXSec(Epion, true, isChargedPion);
86  const double elastic = PionNucleonXSec(Epion, false, isChargedPion);
87  return (total - elastic);
88 }
generate_datataset.total
total
Definition: generate_datataset.py:208
genie::utils::hadxs::berger::PionNucleonXSec
double PionNucleonXSec(double Epion, bool get_total, bool isChargedPion=true)
Definition: HadXSUtils.cxx:96
double genie::utils::hadxs::berger::PionNucleonXSec ( double  Epion,
bool  get_total,
bool  isChargedPion = true 
)

Definition at line 96 of file HadXSUtils.cxx.

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 }
make_norm_csv.output
list output
Definition: make_norm_csv.py:48
genie::constants::kPi0Mass
static const double kPi0Mass
Definition: Constants.h:74
genie::units::mb
static constexpr double mb
Definition: Units.h:79
genie::constants::kPionMass
static const double kPionMass
Definition: Constants.h:73
wirecell.units.pi
float pi
Definition: units.py:11
genie::constants::kPi
static const double kPi
Definition: Constants.h:37
genie::constants::kProtonMass
static const double kProtonMass
Definition: Constants.h:75
genie::constants::kPionMass2
static const double kPionMass2
Definition: Constants.h:86
int genie::utils::hadxs::berger::PionNucleusXSec ( double  tpi,
double  ppistar,
double  t_new,
double  A,
double &  tpilow,
double &  siglow,
double &  tpihigh,
double &  sighigh 
)

Definition at line 205 of file HadXSUtils.cxx.

205  {
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 }
A
#define A
Definition: memgrp.cpp:38
double genie::utils::hadxs::berger::TotalPionNucleonXSec ( double  Epion,
bool  isChargedPion = true 
)

Definition at line 90 of file HadXSUtils.cxx.

92 {
93  return PionNucleonXSec(Epion, true, isChargedPion);
94 }
genie::utils::hadxs::berger::PionNucleonXSec
double PionNucleonXSec(double Epion, bool get_total, bool isChargedPion=true)
Definition: HadXSUtils.cxx:96
Generated by   doxygen 1.8.11