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 91 of file HadXSUtils.cxx.

References PionNucleonXSec(), and produceInfoGainMatrix::total.

Referenced by genie::BergerSehgalFMCOHPiPXSec2015::XSec().

93 {
94  const double total = PionNucleonXSec(Epion, true, isChargedPion);
95  const double elastic = PionNucleonXSec(Epion, false, isChargedPion);
96  return (total - elastic);
97 }
double PionNucleonXSec(double Epion, bool get_total, bool isChargedPion=true)
Definition: HadXSUtils.cxx:105
double genie::utils::hadxs::berger::PionNucleonXSec ( double  Epion,
bool  get_total,
bool  isChargedPion = true 
)

Definition at line 105 of file HadXSUtils.cxx.

References datagram_client::arg1, genie::constants::kPi, genie::constants::kPi0Mass, genie::constants::kPionMass, genie::constants::kPionMass2, genie::constants::kProtonMass, genie::units::mb, confusionMatrixTree::out, runNovaSAM::output, python.hepunit::pi, and ana::Sqrt().

Referenced by InelasticPionNucleonXSec(), TotalPionNucleonXSec(), and genie::BergerSehgalCOHPiPXSec2015::XSec().

107 {
108  // Convert inputs from Genie to those expected by Berger's code:
109 
110  double mpi = kPionMass;
111  double mpi2 = kPionMass2;
112  if (!isChargedPion) {
113  mpi = kPi0Mass;
114  mpi2 = mpi * mpi;
115  }
116 
117  double Epion2 = TMath::Power(Epion,2);
118  double ppi = TMath::Sqrt( TMath::Max(0., Epion2 - mpi2) );
119 
120  if( ppi <= 0.0 ) return 0.0;
121 
122  int out;
123  if( get_total ) out = 0; // Total pion-nucleon cross-section
124  else out = 1; // Elastic pion-nucleon cross-section
125 
126  const double M_pi = mpi; // TODO: used to be kPionMass prior to charge checks
127  const double M_p = kProtonMass; // TODO: should be kNucleonMass ??
128  const double pi = kPi;
129 
130  // Now this is the Berger's code...
131 
132  double afit=1.0;
133  double afit2=1.9;
134  double afit3=0.27;
135  double afit4=0.34;
136  double afit5=0.75;
137  double afit6=1.7;
138 
139  double epi = TMath::Sqrt(M_pi*M_pi + ppi*ppi);
140  double sx = M_p*M_p + M_pi*M_pi + 2.0 * epi * M_p;
141  double Wx = TMath::Sqrt(sx);
142 
143  double s12x = TMath::Sqrt((sx - TMath::Power((M_pi+M_p),2)) * (sx - (M_pi-M_p)*(M_pi-M_p)));
144  double ppistx = s12x/2.0/Wx;
145 
146  double Wdel = 1.232;
147  double arg1 = (Wdel*Wdel - (M_p + M_pi)*(M_p + M_pi))*(Wdel*Wdel - (M_p -M_pi)*(M_p -M_pi));
148  if(arg1<=0) arg1 = 0.0;
149  double gamx = 0.15 * TMath::Power((6.0*ppistx),3) / (1.0 + (6*ppistx)*(6*ppistx));
150  double bwnrx = 1.0 / (4.0*(Wx - Wdel)*(Wx - Wdel) + gamx*gamx);
151  double f1x = afit * 0.3892 * 8.0 * pi / (ppistx*ppistx) * gamx*gamx * bwnrx;
152  double f4x = afit4 * 0.3892 * 8.0 * pi / (ppistx*ppistx) * gamx*gamx * bwnrx;
153 
154  double Wres2 = 1.98;
155  double Gres2 = 0.6;
156  double arg2 = (Wres2*Wres2 - (M_p+M_pi)*(M_p+M_pi)) * (Wres2*Wres2 - (M_p - M_pi)*(M_p - M_pi));
157  if(arg2<=0) arg2 = 0.0;
158  double pp2 = TMath::Sqrt(arg2)/2.0/Wres2;
159  double gam2x = Gres2 * TMath::Power((ppistx/pp2),3);
160  double bwnr2x = 1.0 / (4.0 * (Wx-Wres2)*(Wx-Wres2) + gam2x*gam2x);
161  double f2x = afit2 * 0.3892 * 8.0 * pi / (ppistx*ppistx) * gam2x*gam2x * bwnr2x;
162 
163  double Wres3 = 1.7;
164  double Gres3 = 0.27;
165  double arg3 = (Wres3*Wres3 - (M_p+M_pi)*(M_p+M_pi)) * (Wres3*Wres3 - (M_p - M_pi)*(M_p - M_pi));
166  if(arg3<=0) arg3 = 0.0;
167  double pp3 = TMath::Sqrt(arg3)/2.0/Wres3;
168  double gam3x = Gres3 * TMath::Power((ppistx/pp3),3);
169  double bwnr3x = 1.0 / (4.0 * (Wx-Wres3)*(Wx-Wres3) + gam3x*gam3x);
170  double f3x = afit3 * 0.3892 * 8.0 * pi / (ppistx*ppistx) * gam3x*gam3x * bwnr3x;
171 
172  double Wres5 = 1.52;
173  double Gres5 = 0.10;
174  double arg5 = (Wres5*Wres5 - (M_p+M_pi)*(M_p+M_pi)) * (Wres5*Wres5 - (M_p - M_pi)*(M_p - M_pi));
175  if(arg5<=0) arg5 = 0.0;
176  double pp5 = TMath::Sqrt(arg5)/2.0/Wres5;
177  double gam5x = Gres5 * (ppistx/pp5);
178  double bwnr5x = 1.0 / (4.0 * (Wx-Wres5)*(Wx-Wres5) + gam5x*gam5x);
179  double f5x = afit5 * 0.3892 * 8.0 * pi / (ppistx*ppistx) * gam5x*gam5x * bwnr5x;
180 
181  double Wres6 = 1.69;
182  double Gres6 = 0.140;
183  double arg6 = (Wres6*Wres6 - (M_p+M_pi)*(M_p+M_pi)) * (Wres6*Wres6 - (M_p - M_pi)*(M_p - M_pi));
184  if(arg6<=0) arg6 = 0.0;
185  double pp6 = TMath::Sqrt(arg6)/2.0/Wres6;
186  double gam6x = Gres6 * (ppistx/pp6);
187  double bwnr6x = 1.0 / (4.0 * (Wx-Wres6)*(Wx-Wres6) + gam6x*gam6x);
188  double f6x = afit6 * 0.3892 * 8.0 * pi / (ppistx*ppistx) * gam6x*gam6x * bwnr6x;
189 
190  double output = -9999.99;
191  if(out==0){
192  double bgplusx = 0.834 * ppi + 1.77 * ppi*ppi;
193  double sigplusx = f1x + 0.86 * f2x + 0.7 * f3x + bgplusx;
194  if(ppi>1.9) sigplusx = 20.1 + 15.4 / TMath::Sqrt(ppi);
195  double bgminx = 17.8 * ppi + 6.22 * ppi*ppi - 3.1 * ppi*ppi*ppi;
196  double sigminx = 0.94 * f4x + 0.65 *f5x + 0.61 * f6x + bgminx;
197  if(ppi>1.2) sigminx = 21.64 + 17.29 / TMath::Sqrt(ppi);
198  double sigma_t = (sigplusx + sigminx) / 2.0;
199  output = sigma_t;
200  }else if(out==1){
201  double sgpelx = f1x + 0.411 * f2x;
202  if(ppi>1.9) sgpelx = 2.38 + 17.55 / ppi;
203  double sgmelx = 0.275 * f4x + 0.268 * f5x + 0.281 * f6x + 11.83 * ppi - 3.39 * ppi*ppi;
204  if(ppi>1.5) sgmelx = 3.4 + 11.35 / ppi;
205  double sigma_el = (sgpelx + sgmelx) / 2.0;
206  output = sigma_el;
207  }
208 
209  // Berger's code over, convert to Genie units and return
210  return (output * units::mb);
211 }
const double kPi
const Var kPi0Mass
static const double kPionMass
Definition: Constants.h:74
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:324
static const double kProtonMass
Definition: Constants.h:76
static const double kPionMass2
Definition: Constants.h:87
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 214 of file HadXSUtils.cxx.

Referenced by genie::BergerSehgalFMCOHPiPXSec2015::XSec(), and genie::BergerSehgalCOHPiPXSec2015::XSec().

214  {
215 
216  //Berger code for Tpi<=1.0
217  //Returns the entire dsigma(pi + N -> pi + N)/dt term based on pi-Carbon scattering data
218  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};
219  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};
220  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};
221  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};
222 
223  if(tpi>binedges[12]) return 1;
224  int btu = 1;
225  while(true){
226  if(tpi>binedges[btu]) btu++;
227  else break;
228  }
229  btu--;
230  if(btu<0) btu = 0;
231 
232  tpilow = binedges[btu];
233  tpihigh = binedges[btu + 1];
234 
235  double dsigdzlow = -9999.99;
236  if(btu==0) dsigdzlow = 0.0;
237  else dsigdzlow = 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  siglow = dsigdzlow;
241 
242  double dsigdzhigh = -9999.99;
243  btu++;
244  if(btu>12) btu = 12;
245  if(btu==0) dsigdzhigh = 0.0;
246  else dsigdzhigh = 2.0 * factors[btu]*factors[btu]
247  * parones[btu]*TMath::Power(A/12.0,1.3333333)
248  * TMath::Exp( -1.0 * partwos[btu]*TMath::Power(A/12.0,0.6666666) * t_new * factors[btu]*factors[btu] / (ppistar*ppistar) );
249  sighigh = dsigdzhigh;
250 
251  return 0;
252 }
static const double A
Definition: Units.h:82
double genie::utils::hadxs::berger::TotalPionNucleonXSec ( double  Epion,
bool  isChargedPion = true 
)

Definition at line 99 of file HadXSUtils.cxx.

References PionNucleonXSec().

Referenced by genie::BergerSehgalFMCOHPiPXSec2015::XSec().

101 {
102  return PionNucleonXSec(Epion, true, isChargedPion);
103 }
double PionNucleonXSec(double Epion, bool get_total, bool isChargedPion=true)
Definition: HadXSUtils.cxx:105