41 double Epion2 = TMath::Power(Epion,2);
46 double log10P = TMath::Log10(P);
55 xs =
kInelSig[N-1] + delta * (log10P-log10Pn);
72 double Epion2 = TMath::Power(Epion,2);
77 double log10P = TMath::Log10(P);
86 xs =
kTotSig[N-1] + delta * (log10P-log10Pn);
96 return (total - elastic);
112 if (!isChargedPion) {
117 double Epion2 = TMath::Power(Epion,2);
118 double ppi =
TMath::Sqrt( TMath::Max(0., Epion2 - mpi2) );
120 if( ppi <= 0.0 )
return 0.0;
123 if( get_total ) out = 0;
126 const double M_pi = mpi;
128 const double pi =
kPi;
140 double sx = M_p*M_p + M_pi*M_pi + 2.0 * epi * M_p;
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;
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;
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;
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;
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;
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;
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;
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;
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;
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;
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;
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;
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};
223 if(tpi>binedges[12])
return 1;
226 if(tpi>binedges[btu]) btu++;
232 tpilow = binedges[btu];
233 tpihigh = binedges[btu + 1];
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) );
242 double dsigdzhigh = -9999.99;
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;
static const double kTotdLog10P
static const double kPi0Mass
double TotalPionNucleonXSec(double Epion, bool isChargedPion=true)
double PionNucleonXSec(double Epion, bool get_total, bool isChargedPion=true)
static const double kTotMinLog10P
#define P(a, b, c, d, e, x)
static const double kInelMinLog10P
static const double kIneldLog10P
static const int kInelNDataPoints
double InelasticPionNucleonXSec(double Epion, bool isChargedPion=true)
static const double kTotSig[kTotNDataPoints]
static const double kPionMass
Var Sqrt(const Var &v)
Use to take sqrt of a var.
double InelasticPionNucleonXSec(double Epion, bool isChargedPion=true)
static const double kInelSig[kInelNDataPoints]
static const double kProtonMass
double TotalPionNucleonXSec(double Epion, bool isChargedPion=true)
int PionNucleusXSec(double tpi, double ppistar, double t_new, double A, double &tpilow, double &siglow, double &tpihigh, double &sighigh)
static const double kPionMass2