GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
BostedChristyEMPXSec.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  or see $GENIE/LICENSE
6 
7  For the class documentation see the corresponding header file.
8 */
9 //_________________________________________________________________________
10 #include <vector>
11 #include <string>
12 #include <sstream>
13 
14 #include <TMath.h>
15 #include <TDecayChannel.h>
16 
17 #include "Framework/Algorithm/AlgFactory.h"
18 #include "Framework/Algorithm/AlgConfigPool.h"
19 #include "Framework/Conventions/GBuild.h"
20 #include "Framework/Messenger/Messenger.h"
21 #include "Framework/Conventions/Constants.h"
22 #include "Framework/Utils/KineUtils.h"
23 #include "Framework/ParticleData/PDGLibrary.h"
24 #include "Framework/ParticleData/PDGCodes.h"
25 #include "Framework/ParticleData/PDGUtils.h"
26 #include "Framework/ParticleData/BaryonResUtils.h"
27 #include "Physics/XSectionIntegration/XSecIntegratorI.h"
28 #include "Physics/Resonance/XSection/BostedChristyEMPXSec.h"
29 
30 using namespace genie;
31 using namespace genie::constants;
32 using namespace genie::utils;
33 
34 //____________________________________________________________________________
35 BostedChristyEMPXSec::BostedChristyEMPXSec() :
36 XSecAlgorithmI("genie::BostedChristyEMPXSec")
37 {
38 }
39 //____________________________________________________________________________
40 BostedChristyEMPXSec::BostedChristyEMPXSec(string config) :
41 XSecAlgorithmI("genie::BostedChristyEMPXSec", config)
42 {
43 }
44 //____________________________________________________________________________
45 BostedChristyEMPXSec::~BostedChristyEMPXSec()
46 {
47 
48 }
49 //____________________________________________________________________________
50 // resonance cross section fit of electron-proton and electron-deuterium scattering data
51 // sf=0 \sigma_T
52 // sf=1 \sigma_L
53 double BostedChristyEMPXSec::sigmaR(int sf, double Q2, double W, bool isDeuterium=false) const
54 {
55  const std::array<std::array<double, 3>, 7> &fBR = !isDeuterium?fBRp:fBRD;
56  const std::array<std::array<double, 4>, 7> &fRescoefT = !isDeuterium?fRescoefTp:fRescoefTD;
57  if (isDeuterium)
58  sf=0;
59 
60  double W2 = W*W;
61  // proton mass
62  double Mp = fMP;
63  // pion mass
64  double Mpi = fMpi0;
65  // eta-meson mass
66  double Meta = fMeta;
67 
68  double Mp2 = Mp*Mp;
69  double Mpi2 = Mpi*Mpi;
70  double Meta2 = Meta*Meta;
71 
72  // Calculate kinematics needed for threshold Relativistic B-W
73 
74  // Ref.1, Eq. (10)
75  double k = (W2 - Mp2)/2./Mp;
76  // Ref.1, Eq. (11)
77  double kcm = (W2 - Mp2)/2./W;
78  // mesons energy and momentim
79  double Epicm = (W2 + Mpi2 - Mp2)/2./W; // pion energy in CMS
80  double ppicm = TMath::Sqrt(TMath::Max(0.0, Epicm*Epicm - Mpi2)); // pion momentum in CMS
81  double Epi2cm = (W2 + 4*Mpi2 - Mp2)/2./W; // two pion energy in CMS
82  double ppi2cm = TMath::Sqrt(TMath::Max(0.0, Epi2cm*Epi2cm - 4*Mpi2)); // two pion energi n CMS
83  double Eetacm = (W2 + Meta2 - Mp2 )/2./W; // eta energy in CMS
84  double petacm = TMath::Sqrt(TMath::Max(0.0, Eetacm*Eetacm - Meta2)); // eta energy in CMS
85 
86  double xsec = 0.0;
87 
88  // going through seven resonances
89  for (int i=0;i<7;i++)
90  {
91  // resonance mass squared
92  double MassRes2 = fMassRes[i]*fMassRes[i];
93  // resonance damping parameter
94  double x0 = i!=0?0.215:!isDeuterium?0.14462:0.1446;
95  // Ref.1, Eq. (12)
96  double kr = (MassRes2-Mp2)/2./Mp;
97  // Ref.1, Eq. (13)
98  double kcmr = (MassRes2-Mp2)/2./fMassRes[i];
99 
100  // formulas analogous to the above with substitution W->MR_i
101  double Epicmr = (MassRes2 + Mpi2 - Mp2)/2./fMassRes[i];
102  double ppicmr = TMath::Sqrt(TMath::Max(0.0, Epicmr*Epicmr - Mpi2));
103  double Epi2cmr = (MassRes2 + 4.*Mpi2 - Mp2)/2./fMassRes[i];
104  double ppi2cmr = TMath::Sqrt(TMath::Max(0.0, Epi2cmr*Epi2cmr - 4.*Mpi2));
105  double Eetacmr = (MassRes2 + Meta2 - Mp2)/2./fMassRes[i];
106  double petacmr = TMath::Sqrt(TMath::Max(0.0, Eetacmr*Eetacmr - Meta2));
107 
108  // Calculate partial widths
109  // Ref.1, Eq. (15) for single pion
110  double pwid0 = fWidthRes[i]*TMath::Power(ppicm/ppicmr, 1.+2.*fAngRes[i])*
111  TMath::Power((ppicmr*ppicmr + x0*x0)/(ppicm*ppicm+x0*x0), fAngRes[i]); // 1-pion decay mode
112  // Ref.1, Eq. (16) for two pions
113  double pwid1 = 0;
114  if (!isDeuterium || (isDeuterium && i!=1))
115  pwid1 = W/fMassRes[i]*fWidthRes[i]*TMath::Power(ppi2cm/ppi2cmr, 4.+2.*fAngRes[i])*
116  TMath::Power((ppi2cmr*ppi2cmr + x0*x0)/(ppi2cm*ppi2cm + x0*x0), 2.+fAngRes[i]); // 2-pion decay mode
117  else
118  pwid1 = fWidthRes[i]*TMath::Power(petacm/petacmr, 1.+2.*fAngRes[i])*TMath::Power((ppi2cmr*ppi2cmr + x0*x0)/(ppi2cm*ppi2cm + x0*x0), fAngRes[i]);
119 
120 
121  double pwid2 = 0.; // eta decay mode
122  // Ref.1, Eq. (15) for eta
123  if(!isDeuterium && (i==1 || i==4))
124  pwid2 = fWidthRes[i]*TMath::Power(petacm/petacmr, 1.+2.*fAngRes[i])*TMath::Power((petacmr*petacmr + x0*x0)/(petacm*petacm + x0*x0), fAngRes[i]); // eta decay only for S11's
125 
126  // Ref.1, Eq. (17)
127  double pgam = fWidthRes[i]*(kcm/kcmr)*(kcm/kcmr)*(kcmr*kcmr+x0*x0)/(kcm*kcm+x0*x0);
128  // Ref.1, Eq. (14)
129  double width = fBR[i][0]*pwid0+fBR[i][1]*pwid1+fBR[i][2]*pwid2;
130 
131  // Begin resonance Q^2 dependence calculations
132  double A;
133 
134  if (sf==0)
135  // Ref.1, Eq. (18)
136  A = fRescoefT[i][0]*(1.+fRescoefT[i][1]*Q2/(1.+fRescoefT[i][2]*Q2))/TMath::Power(1.+Q2/0.91, fRescoefT[i][3]);
137  else
138  // Ref.1, Eq. (19)
139  A = fRescoefL[i][0]*Q2/(1.+fRescoefL[i][1]*Q2)*TMath::Exp(-1.*fRescoefL[i][2]*Q2);
140 
141 
142  // Ref.1, Eq. (9)
143  double BW = kr/k*kcmr/kcm/fWidthRes[i]*width*pgam/((W2 - MassRes2)*(W2 - MassRes2) + MassRes2*width*width);
144 
145  // Ref.1, Eq. (8) divided by W
146  xsec += BW*A*A;
147  }
148  // restore factor W in Ref.1, Eq. (8)
149  return W*xsec*units::ub;
150 }
151 //____________________________________________________________________________
152 // nonresonance cross section fit of electron-proton and electron-deuterium scattering data
153 // sf=0 \sigma_T
154 // sf=1 \sigma_L
155 double BostedChristyEMPXSec::sigmaNR(int sf, double Q2, double W, bool isDeuterium=false) const
156 {
157  const std::array<std::array<double, 5>, 2> &fNRcoefT = !isDeuterium?fNRcoefTp:fNRcoefTD;
158  if (isDeuterium)
159  sf=0;
160  double W2 = W*W;
161  double Mp = fMP;
162  double Mpi = fMpi0;
163  double Mp2 = Mp*Mp;
164 
165  double Wdif = W - (Mp + Mpi);
166 
167  double m0 = (sf==0) ? 0.125 : 4.2802; //Ref.1, Eqs.(22, 24)
168 
169  double Q20 = (sf==0) ? 0.05 : 0.125; //Ref.1, Eqs.(22, 24)
170 
171  double xpr = 1./(1.+(W2-(Mp+Mpi)*(Mp+Mpi))/(Q2+Q20)); // Ref.1, Eq.(22)
172 
173  double xsec = 0.0;
174 
175 
176  if (sf==0)
177  {
178 
179  for (int i=0;i<2;i++)
180  {
181  double h_nr = fNRcoefT[i][0]/TMath::Power(Q2+fNRcoefT[i][1], fNRcoefT[i][2]+fNRcoefT[i][3]*Q2+fNRcoefT[i][4]*Q2*Q2); // Ref.1, Eq. (21)
182  xsec += h_nr*TMath::Power(Wdif, 1.5+i);
183  }
184 
185  xsec *= xpr;
186  }
187  else
188  {
189  double xb = Q2/(Q2+W2-Mp2);
190  double t = TMath::Log(TMath::Log((Q2+m0)/0.330/0.330)/TMath::Log(m0/0.330/0.330)); // Ref.1, Eq. (24)
191  xsec += fNRcoefL[0]*TMath::Power(1.-xpr, fNRcoefL[2]+fNRcoefL[1]*t)/(1.-xb)/(Q2+Q20)
192  *TMath::Power(Q2/(Q2+Q20), fNRcoefL[3])*TMath::Power(xpr, fNRcoefL[4]+fNRcoefL[5]*t); // Ref.1, Eq. (23)
193  }
194 
195  return xsec*units::ub;
196 }
197 //___________________________________________________________________________
198 // Calculate proton and neutron with Fermi smearing of a nulei
199 void BostedChristyEMPXSec::FermiSmearingA(double Q2, double W, double pF, double Es, double & F1p, double & F1d, double & sigmaT, double & sigmaL) const
200 {
201  // The numbers in arrays bellow were not supposed to change in the original
202  // fortran code and therefore are not configurable
203  static constexpr std::array<double, 99> fyp
204  {0.0272,0.0326,0.0390,0.0464,0.0551,0.0651,0.0766,0.0898,0.1049,
205  0.1221,0.1416,0.1636,0.1883,0.2159,0.2466,0.2807,0.3182,0.3595,
206  0.4045,0.4535,0.5066,0.5637,0.6249,0.6901,0.7593,0.8324,0.9090,
207  0.9890,1.0720,1.1577,1.2454,1.3349,1.4254,1.5163,1.6070,1.6968,
208  1.7849,1.8705,1.9529,2.0313,2.1049,2.1731,2.2350,2.2901,2.3379,
209  2.3776,2.4090,2.4317,2.4454,2.4500,2.4454,2.4317,2.4090,2.3776,
210  2.3379,2.2901,2.2350,2.1731,2.1049,2.0313,1.9529,1.8705,1.7849,
211  1.6968,1.6070,1.5163,1.4254,1.3349,1.2454,1.1577,1.0720,0.9890,
212  0.9090,0.8324,0.7593,0.6901,0.6249,0.5637,0.5066,0.4535,0.4045,
213  0.3595,0.3182,0.2807,0.2466,0.2159,0.1883,0.1636,0.1416,0.1221,
214  0.1049,0.0898,0.0766,0.0651,0.0551,0.0464,0.0390,0.0326,0.0272};
215 
216  static constexpr std::array<double, 99> xxp
217  {-3.000,-2.939,-2.878,-2.816,-2.755,-2.694,-2.633,-2.571,-2.510,
218  -2.449,-2.388,-2.327,-2.265,-2.204,-2.143,-2.082,-2.020,-1.959,
219  -1.898,-1.837,-1.776,-1.714,-1.653,-1.592,-1.531,-1.469,-1.408,
220  -1.347,-1.286,-1.224,-1.163,-1.102,-1.041,-0.980,-0.918,-0.857,
221  -0.796,-0.735,-0.673,-0.612,-0.551,-0.490,-0.429,-0.367,-0.306,
222  -0.245,-0.184,-0.122,-0.061, 0.000, 0.061, 0.122, 0.184, 0.245,
223  0.306, 0.367, 0.429, 0.490, 0.551, 0.612, 0.673, 0.735, 0.796,
224  0.857, 0.918, 0.980, 1.041, 1.102, 1.163, 1.224, 1.286, 1.347,
225  1.408, 1.469, 1.531, 1.592, 1.653, 1.714, 1.776, 1.837, 1.898,
226  1.959, 2.020, 2.082, 2.143, 2.204, 2.265, 2.327, 2.388, 2.449,
227  2.510, 2.571, 2.633, 2.694, 2.755, 2.816, 2.878, 2.939, 3.000};
228 
229  double MN = fPM;
230  double MN2 = MN*MN;
231  double Mp = fMP;
232  double Mp2 = Mp*Mp;
233  double W2 = W*W;
234 
235  double nu = (W2 - MN2 + Q2)/2./MN;
236  double qv = TMath::Sqrt(nu*nu + Q2);
237  // assume this is 2*pf*qv
238  double dW2dpF = 2.*qv;
239  double dW2dEs = 2.*(nu + MN);
240  // switched to using 99 bins!
241  F1p = 0;
242  F1d = 0;
243  sigmaT = 0;
244  sigmaL = 0;
245  for (int i=0; i<99; i++)
246  {
247  double fyuse = fyp[i]/100.;
248  double W2p = W2 + xxp[i]*pF*dW2dpF - Es*dW2dEs;
249  if(W2p>1.159)
250  {
251  //proton
252  double Wp = TMath::Sqrt(W2p);
253  double sigmaTp = sigmaR(0, Q2, Wp) + sigmaNR(0, Q2, Wp);
254  double sigmaLp = sigmaR(1, Q2, Wp) + sigmaNR(1, Q2, Wp);
255  double F1pp = sigmaTp*(W2p-Mp2)/8./kPi2/kAem;
256  //neutron
257  double sigmaTd = sigmaR(0, Q2, Wp, true) + sigmaNR(0, Q2, Wp, true);
258  double F1dp = sigmaTd*(W2p-Mp2)/8./kPi2/kAem;
259  F1d += F1dp*fyuse;
260  F1p += F1pp*fyuse;
261  sigmaT += sigmaTp*fyuse;
262  sigmaL += sigmaLp*fyuse;
263  }
264  }
265 
266 }
267 //___________________________________________________________________________
268 // Calculate proton and neutron with Fermi smearing of a deuteron
269 void BostedChristyEMPXSec::FermiSmearingD(double Q2, double W, double & F1, double & R, double & sigmaT, double & sigmaL, bool isDeuterium=false) const
270 {
271  // The numbers in arrays bellow were not supposed to change in the original
272  // fortran code and therefore are not configurable
273  static constexpr std::array<double, 20> fyd
274  {0.4965, 0.4988, 0.4958, 0.5008, 0.5027, 0.5041, 0.5029, 0.5034,
275  0.4993, 0.5147, 0.5140, 0.4975, 0.5007, 0.4992, 0.4994, 0.4977,
276  0.5023, 0.4964, 0.4966, 0.4767};
277 
278  static constexpr std::array<double, 20> avpz
279  {-0.1820,-0.0829,-0.0590,-0.0448,-0.0345,-0.0264,-0.0195, -0.0135,
280  -0.0079,-0.0025, 0.0029, 0.0083, 0.0139, 0.0199, 0.0268, 0.0349,
281  0.0453, 0.0598, 0.0844, 0.1853};
282 
283  static constexpr std::array<double, 20> avp2
284  {0.0938, 0.0219, 0.0137, 0.0101, 0.0081, 0.0068, 0.0060, 0.0054,
285  0.0051, 0.0049, 0.0050, 0.0051, 0.0055, 0.0060, 0.0069, 0.0081,
286  0.0102, 0.0140, 0.0225, 0.0964};
287 
288  // Look up tables for deuteron in fine bins for sub threshold
289  static constexpr std::array<double, 200> fydf
290  {0.00001,0.00002,0.00003,0.00005,0.00006,0.00009,0.00010,0.00013,
291  0.00015,0.00019,0.00021,0.00026,0.00029,0.00034,0.00038,0.00044,
292  0.00049,0.00057,0.00062,0.00071,0.00078,0.00089,0.00097,0.00109,
293  0.00119,0.00134,0.00146,0.00161,0.00176,0.00195,0.00211,0.00232,
294  0.00252,0.00276,0.00299,0.00326,0.00352,0.00383,0.00412,0.00447,
295  0.00482,0.00521,0.00560,0.00603,0.00648,0.00698,0.00747,0.00803,
296  0.00859,0.00921,0.00985,0.01056,0.01126,0.01205,0.01286,0.01376,
297  0.01467,0.01569,0.01671,0.01793,0.01912,0.02049,0.02196,0.02356,
298  0.02525,0.02723,0.02939,0.03179,0.03453,0.03764,0.04116,0.04533,
299  0.05004,0.05565,0.06232,0.07015,0.07965,0.09093,0.10486,0.12185,
300  0.14268,0.16860,0.20074,0.24129,0.29201,0.35713,0.44012,0.54757,
301  0.68665,0.86965,1.11199,1.43242,1.86532,2.44703,3.22681,4.24972,
302  5.54382,7.04016,8.48123,9.40627,9.40627,8.48123,7.04016,5.54382,
303  4.24972,3.22681,2.44703,1.86532,1.43242,1.11199,0.86965,0.68665,
304  0.54757,0.44012,0.35713,0.29201,0.24129,0.20074,0.16860,0.14268,
305  0.12185,0.10486,0.09093,0.07965,0.07015,0.06232,0.05565,0.05004,
306  0.04533,0.04116,0.03764,0.03453,0.03179,0.02939,0.02723,0.02525,
307  0.02356,0.02196,0.02049,0.01912,0.01793,0.01671,0.01569,0.01467,
308  0.01376,0.01286,0.01205,0.01126,0.01056,0.00985,0.00921,0.00859,
309  0.00803,0.00747,0.00698,0.00648,0.00603,0.00560,0.00521,0.00482,
310  0.00447,0.00412,0.00383,0.00352,0.00326,0.00299,0.00276,0.00252,
311  0.00232,0.00211,0.00195,0.00176,0.00161,0.00146,0.00134,0.00119,
312  0.00109,0.00097,0.00089,0.00078,0.00071,0.00062,0.00057,0.00049,
313  0.00044,0.00038,0.00034,0.00029,0.00026,0.00021,0.00019,0.00015,
314  0.00013,0.00010,0.00009,0.00006,0.00005,0.00003,0.00002,0.00001};
315 
316  static constexpr std::array<double, 200> avp2f
317  {1.0,0.98974,0.96975,0.96768,0.94782,0.94450,0.92494,0.92047,
318  0.90090,0.89563,0.87644,0.87018,0.85145,0.84434,0.82593,0.81841,
319  0.80021,0.79212,0.77444,0.76553,0.74866,0.73945,0.72264,0.71343,
320  0.69703,0.68740,0.67149,0.66182,0.64631,0.63630,0.62125,0.61154,
321  0.59671,0.58686,0.57241,0.56283,0.54866,0.53889,0.52528,0.51581,
322  0.50236,0.49291,0.47997,0.47063,0.45803,0.44867,0.43665,0.42744,
323  0.41554,0.40656,0.39511,0.38589,0.37488,0.36611,0.35516,0.34647,
324  0.33571,0.32704,0.31656,0.30783,0.29741,0.28870,0.27820,0.26945,
325  0.25898,0.25010,0.23945,0.23023,0.21943,0.20999,0.19891,0.18911,
326  0.17795,0.16793,0.15669,0.14667,0.13553,0.12569,0.11504,0.10550,
327  0.09557,0.08674,0.07774,0.06974,0.06184,0.05484,0.04802,0.04203,
328  0.03629,0.03129,0.02654,0.02247,0.01867,0.01545,0.01251,0.01015,
329  0.00810,0.00664,0.00541,0.00512,0.00512,0.00541,0.00664,0.00810,
330  0.01015,0.01251,0.01545,0.01867,0.02247,0.02654,0.03129,0.03629,
331  0.04203,0.04802,0.05484,0.06184,0.06974,0.07774,0.08674,0.09557,
332  0.10550,0.11504,0.12569,0.13553,0.14667,0.15669,0.16793,0.17795,
333  0.18911,0.19891,0.20999,0.21943,0.23023,0.23945,0.25010,0.25898,
334  0.26945,0.27820,0.28870,0.29741,0.30783,0.31656,0.32704,0.33571,
335  0.34647,0.35516,0.36611,0.37488,0.38589,0.39511,0.40656,0.41554,
336  0.42744,0.43665,0.44867,0.45803,0.47063,0.47997,0.49291,0.50236,
337  0.51581,0.52528,0.53889,0.54866,0.56283,0.57241,0.58686,0.59671,
338  0.61154,0.62125,0.63630,0.64631,0.66182,0.67149,0.68740,0.69703,
339  0.71343,0.72264,0.73945,0.74866,0.76553,0.77444,0.79212,0.80021,
340  0.81841,0.82593,0.84434,0.85145,0.87018,0.87644,0.89563,0.90090,
341  0.92047,0.92494,0.94450,0.94782,0.96768,0.96975,0.98974,1.0};
342 
343  double W2=W*W;
344  double MN = fAM;
345  double MN2 = MN*MN;
346  double MD = fMD;
347  double Mp = fMP;
348  double Mp2 = Mp*Mp;
349  double nu = (W2 - MN2 + Q2)/2./MN;
350  double qv = TMath::Sqrt(nu*nu + Q2);
351  F1 = 0.;
352  R = 0.;
353  sigmaT = 0.;
354  sigmaL = 0.;
355  // Do fast 20 bins if abvoe threshold
356  if(W2>1.30)
357  {
358  for (int ism = 0; ism<20; ism++)
359  {
360  double W2p = TMath::Power(MD + nu - sqrt(MN2 + avp2[ism]),2) - qv*qv + 2.*qv*avpz[ism] - avp2[ism];
361  if(W2p>1.155)
362  {
363  double Wp = TMath::Sqrt(W2p);
364  double sigtp = sigmaR(0, Q2, Wp, isDeuterium) + sigmaNR(0, Q2, Wp, isDeuterium);
365  double F1p = sigtp*(W2p-Mp2)/8./kPi2/kAem;
366  F1 += F1p*fyd[ism]/10.;
367  if (!isDeuterium)
368  {
369  double siglp = sigmaR(1, Q2, Wp) + sigmaNR(1, Q2, Wp);
370  sigmaL += siglp*fyd[ism]/10.;
371  sigmaT += sigtp*fyd[ism]/10.;
372  }
373  }
374  }
375  }
376  else
377  {
378  for (int ism = 0;ism<200;ism++)
379  {
380  double pz = -1. + 0.01*(ism + 0.5);
381  // Need avp2f term to get right behavior x>1!
382  double W2p = TMath::Power(MD + nu - sqrt(MN2 + avp2f[ism]),2) - qv*qv + 2.*qv*pz - avp2f[ism];
383  if(W2p>1.155)
384  {
385  double Wp = TMath::Sqrt(W2p);
386  double sigtp = sigmaR(0, Q2, Wp, isDeuterium) + sigmaNR(0, Q2, Wp, isDeuterium);
387  double F1p = sigtp*(W2p-Mp2)/8./kPi2/kAem;
388  F1 += F1p*fydf[ism]/100.;
389  if (!isDeuterium)
390  {
391  double siglp = sigmaR(1, Q2, Wp) + sigmaNR(1, Q2, Wp);
392  sigmaT += sigtp*fydf[ism]/100.;
393  sigmaL += siglp*fydf[ism]/100.;
394  }
395  }
396  }
397  }
398  if (isDeuterium && fUseMEC)
399  // Ref.2, Eq. (20)
400  F1 += fMECcoef[0]*TMath::Exp(-(W - fMECcoef[1])*(W - fMECcoef[1])/fMECcoef[2])/
401  TMath::Power(1. + TMath::Max(0.3,Q2)/fMECcoef[3],fMECcoef[4])*TMath::Power(nu, fMECcoef[5]);
402  if(!isDeuterium && sigmaT!=0.)
403  R = sigmaL/sigmaT;
404 
405 }
406 //____________________________________________________________________________
407 double BostedChristyEMPXSec::MEC2009(int A, double Q2, double W) const
408 {
409  double F1 = 0.0;
410  double W2 = W*W;
411  double Mp = fAM;
412  double Mp2 = Mp*Mp;
413  if(W2<=0.0)
414  return F1;
415  double nu = (W2 - Mp2 + Q2)/2./Mp;
416  double x = Q2/2.0/Mp/nu;
417 
418  if(A<=2)
419  return F1;
420 
421  double p18;
422  for (const auto& kv : fMEC2009p18)
423  {
424  p18 = kv.second;
425  if (A<=kv.first)
426  break;
427  }
428 
429  F1 = fMEC2009coef[0]*TMath::Exp(-(W - fMEC2009coef[1])*(W - fMEC2009coef[1])/fMEC2009coef[2])/
430  TMath::Power(1. + TMath::Max(0.3,Q2)/fMEC2009coef[3],fMEC2009coef[4])*TMath::Power(nu, fMEC2009coef[5])*(1.0 +
431  p18*TMath::Power(A, fMEC2009coef[6] + fMEC2009coef[7]*x));
432 
433  if(F1<=1.0E-9)
434  F1 = 0.0;
435 
436  return F1;
437 }
438 //____________________________________________________________________________
439 double BostedChristyEMPXSec::XSec(
440  const Interaction * interaction, KinePhaseSpace_t kps) const
441 {
442  if(! this -> ValidProcess (interaction) ) return 0.;
443  if(! this -> ValidKinematics (interaction) ) return 0.;
444 
445  // Get kinematical parameters
446  const Kinematics & kinematics = interaction -> Kine();
447  const InitialState & init_state = interaction -> InitState();
448  const Target & target = init_state.Tgt();
449  int A = target.A();
450  int Z = target.Z();
451  double E = init_state.ProbeE(kRfHitNucRest);
452  double W = kinematics.W();
453  double Q2 = kinematics.Q2();
454  double Wsq = W*W;
455  // Cross section for proton or neutron
456 
457  double Mp = fMP;
458  double Mp2 = Mp*Mp;
459  double MN = fPM;
460  double MN2 = MN*MN;
461 
462  double nu = (Wsq - MN2 + Q2)/2./MN;
463  double x = Q2/2./MN/nu;
464 
465  double sigmaT, sigmaL, F1p, R, W1;
466  // Cross section for proton or neutron
467  if (A<2 && Wsq>1.155)
468  {
469  double xb = Q2/(Wsq+Q2-Mp2);
470  sigmaT = sigmaR(0, Q2, W) + sigmaNR(0, Q2, W);
471  sigmaL = sigmaR(1, Q2, W) + sigmaNR(1, Q2, W);
472  F1p = sigmaT*(Wsq-Mp2)/8./kPi2/kAem;
473  R = sigmaL/sigmaT;
474  // If neutron, subtract proton from deuteron. Factor of two to
475  // convert from per nucleon to per deuteron
476  if(Z==0)
477  {
478  sigmaT = sigmaR(0, Q2, W, true) + sigmaNR(0, Q2, W, true);
479  double F1d = sigmaT*(Wsq-Mp2)/8./kPi2/kAem;
480  F1p = 2.*F1d - F1p;
481  }
482  W1 = F1p/MN;
483  }
484 
485  // For deuteron
486  if(A==2)
487  {
488  double Rd, F1c, F1d;
489  //get Fermi-smeared R from Erics proton fit
490  FermiSmearingD(Q2, W, F1c, R, sigmaT, sigmaL);
491  //get fit to F1 in deuteron, per nucleon
492  FermiSmearingD(Q2, W, F1d, Rd, sigmaT, sigmaL, true);
493  //convert to W1 per deuteron
494  W1 = F1d/MN*2.0;
495  }
496 
497  //For nuclei
498  if (A>2)
499  {
500  // Modifed to use Superscaling from Ref. 3
501  double Es, pF, kF;
502  for (const auto& kv : fNucRmvE)
503  {
504  Es = kv.second;
505  if (A<=kv.first)
506  break;
507  }
508  for (const auto& kv : fKFTable)
509  {
510  kF = kv.second;
511  if (A<=kv.first)
512  break;
513  }
514  // adjust pf to give right width based on kf
515  pF = 0.5*kF;
516  double F1d;
517  FermiSmearingA(Q2, W, pF, Es, F1p, F1d, sigmaT, sigmaL);
518  R = 0.;
519  if(sigmaT>0.)
520  R = sigmaL/sigmaT;
521  W1 = (2.*Z*F1d + (A - 2.*Z)*(2.*F1d - F1p))/MN;
522 
523  W1 *= (fAfitcoef[0] + x*(fAfitcoef[1] + x*(fAfitcoef[2] + x*(fAfitcoef[3] + x*(fAfitcoef[4] + x*fAfitcoef[5])))));
524 
525  if(W>0.)
526  W1 *= TMath::Power(fAfitcoef[6] + (fAfitcoef[7]*W + fAfitcoef[8]*Wsq)/(fAfitcoef[9] + fAfitcoef[10]*Q2),2);
527 
528  double F1M = MEC2009(A, Q2, W);
529 
530  W1 += F1M;
531  if(Wsq>0.)
532  R *= (fAfitcoef[11] + fAfitcoef[12]*A);
533  }
534 
535  double emcfac = FitEMC(x, A);
536 
537  W1 *= emcfac;
538 
539  double nu2 = nu*nu;
540  double Eprime = E - nu;
541  double sin2theta_2 = Q2/4/E/Eprime;
542  double cos2theta_2 = 1 - sin2theta_2;
543  double W2 = W1*(1 + R)/ (1+nu2/Q2);
544  double xsec = 4*Eprime*Eprime*kAem2/Q2/Q2*(2*W1*sin2theta_2 + W2*cos2theta_2); // d2xsec/dOmegadEprime
545  double jacobian = W*kPi/E/Eprime/MN;
546  xsec*= jacobian; // d2xsec/dOmegadEprime-> d2xsec/dWdQ2
547 
548  // The algorithm computes d^2xsec/dWdQ2
549  // Check whether variable tranformation is needed
550  if(kps!=kPSWQ2fE) {
551  double J = utils::kinematics::Jacobian(interaction,kPSWQ2fE,kps);
552  xsec *= J;
553  }
554 
555  return xsec;
556 
557 }
558 //____________________________________________________________________________
559 // Fit to EMC effect. Steve Rock 8/3/94
560 // Funciton returns value of sigma(A)/sigma(d)
561 // with no isoscalerity correction
562 // A = atomic number
563 // x = Bjorken x.
564 //
565 // Fit of sigma(A)/sigma(d) to form C*A**alpha where A is atomic number
566 // First data at each x was fit to form C*A**alpha. The point A=2(d)
567 // was includded with a value of 1 and an error of about 2%.
568 // For x>=.125 Javier Gomez fit of 7/93 to E139 data was used.
569 // For .09 >=x>=.0085 NMC data from Amaudruz et al Z. Phys C. 51,387(91)
570 // Steve did the fit for alpha and C to the He/d. C/d and Ca/d NMC data.
571 // Alpha(x) was fit to a 9 term polynomial a0 +a1*x +a2*x**2 + a3*x**3 ..
572 // C(x) was fit to a 3 term polynomial in natural logs as
573 // Ln(C) = c0 + c1*Ln(x) + c2*[Ln(x)]**2.
574 //
575 // 6/2/98 ***** Bug (which set x= .00885 if x was out of range) fixed
576 // also gave value at x=.0085 if x>.88
577 // 11/05 PYB PEB modified to use value at x=0.7 if x>0.7, because
578 // beyond that assume Fermi motion is giving the rise, and we
579 // already are taking that into account with the y-smearing of
580 // the inelastic
581 //____________________________________________________________________________
582 double BostedChristyEMPXSec::FitEMC(double x, int A) const
583 {
584  double fitemc = 1.;
585  if(A<=2)
586  return fitemc;
587 
588  double x_u;
589  if (x>0.70 || x<0.0085)
590  //Out of range of fit
591  {
592  if(x<0.0085)
593  x_u = .0085;
594  if(x>0.70)
595  x_u = 0.70;
596  }
597  else
598  x_u = x;
599 
600  double ln_c = fEMCc[0];
601  for (int i = 1; i<=2; i++)
602  ln_c += fEMCc[i]*TMath::Power(TMath::Log(x_u), i);
603  double c = TMath::Exp(ln_c);
604 
605  double alpha = fEMCalpha[0];
606  for (int i = 1; i<=8; i++)
607  alpha += fEMCalpha[i]*TMath::Power(x_u, i);
608 
609  fitemc = c*TMath::Power(A, alpha);
610  return fitemc;
611 }
612 //____________________________________________________________________________
613 double BostedChristyEMPXSec::Integral(const Interaction * interaction) const
614 {
615  double xsec = fXSecIntegrator->Integrate(this,interaction);
616  return xsec;
617 }
618 //____________________________________________________________________________
619 bool BostedChristyEMPXSec::ValidProcess(const Interaction * interaction) const
620 {
621  if(interaction->TestBit(kISkipProcessChk)) return true;
622 
623  const InitialState & init_state = interaction->InitState();
624  const ProcessInfo & proc_info = interaction->ProcInfo();
625 
626  if(!proc_info.IsResonant() ) return false;
627 
628 
629  int hitnuc = init_state.Tgt().HitNucPdg();
630  bool is_pn = (pdg::IsProton(hitnuc) || pdg::IsNeutron(hitnuc));
631 
632  if (!is_pn) return false;
633 
634  int probe = init_state.ProbePdg();
635  bool is_em = proc_info.IsEM();
636 
637  bool l_em = (pdg::IsChargedLepton(probe) && is_em );
638 
639  if (!l_em) return false;
640 
641  return true;
642 }
643 //____________________________________________________________________________
644 bool BostedChristyEMPXSec::ValidKinematics(const Interaction * interaction) const
645 {
646  const Kinematics & kinematics = interaction -> Kine();
647  double W = kinematics.W();
648  double Q2 = kinematics.Q2();
649  if (W<fWmin || W>fWmax)
650  return false;
651  if (Q2<fQ2min || Q2>fQ2max)
652  return false;
653 
654  return true;
655 }
656 //____________________________________________________________________________
657 void BostedChristyEMPXSec::Configure(const Registry & config)
658 {
659  Algorithm::Configure(config);
660  this->LoadConfig();
661 }
662 //____________________________________________________________________________
663 void BostedChristyEMPXSec::Configure(string config)
664 {
665  Algorithm::Configure(config);
666  this->LoadConfig();
667 }
668 //____________________________________________________________________________
669 void BostedChristyEMPXSec::BranchingRatios(int respdg, double & brpi, double & breta) const
670 {
671  brpi = 0.;
672  breta = 0.;
673  PDGLibrary * pdglib = PDGLibrary::Instance();
674  TParticlePDG * res_pdg = pdglib->Find(respdg);
675  if (res_pdg != 0)
676  {
677  for (int nch = 0; nch < res_pdg->NDecayChannels(); nch++)
678  {
679  TDecayChannel * ch = res_pdg->DecayChannel(nch);
680  if (ch->NDaughters() == 2)
681  {
682  int first_daughter_pdg = ch->DaughterPdgCode (0);
683  int second_daughter_pdg = ch->DaughterPdgCode (1);
684  if ((genie::pdg::IsNucleon(first_daughter_pdg ) && genie::pdg::IsPion(second_daughter_pdg)) ||
685  (genie::pdg::IsNucleon(second_daughter_pdg) && genie::pdg::IsPion(first_daughter_pdg )))
686  brpi += ch->BranchingRatio();
687  if (first_daughter_pdg == kPdgEta || second_daughter_pdg == kPdgEta)
688  breta += ch->BranchingRatio();
689  }
690  }
691  }
692 }
693 //____________________________________________________________________________
694 void BostedChristyEMPXSec::LoadConfig(void)
695 {
696 
697  PDGLibrary * pdglib = PDGLibrary::Instance();
698  GetParamDef("BostedChristyFitEM-PM", fPM, pdglib->Find(kPdgProton)->Mass());
699  GetParamDef("BostedChristyFitEM-MP", fMP, pdglib->Find(kPdgProton)->Mass());
700  GetParamDef("BostedChristyFitEM-AM", fAM, pdglib->Find(kPdgProton)->Mass());
701  GetParamDef("BostedChristyFitEM-MD", fMD, pdglib->Find(kPdgTgtDeuterium)->Mass());
702  GetParamDef("BostedChristyFitEM-Mpi0", fMpi0, pdglib->Find(kPdgPi0)->Mass());
703  GetParamDef("BostedChristyFitEM-Meta", fMeta, pdglib->Find(kPdgEta)->Mass());
704  GetParamDef("BostedChristyFitEM-Wmin", fWmin, 0.0);
705  GetParamDef("BostedChristyFitEM-Wmax", fWmax, 3.0);
706  GetParamDef("BostedChristyFitEM-Q2min", fQ2min, 0.0);
707  GetParamDef("BostedChristyFitEM-Q2max", fQ2max, 10.0);
708  GetParamDef("BostedChristyFitEM-UseMEC", fUseMEC, true);
709 
710  double BRpi, BReta;
711  double brpi, breta;
712 
713  std::vector<double> vBRpi1;
714  std::vector<double> vBRpi2;
715  std::vector<double> vBReta;
716 
717  // load braching ratios for pi
718  bool useBRpi1Default = (GetParamVect("BostedChristyFitEM-PionBRp", vBRpi1, false)<7);
719  bool useBRpi2Default = (GetParamVect("BostedChristyFitEM-PionBRD", vBRpi2, false)<7);
720  // load braching ratios for eta
721  bool useBRetaDefault = (GetParamVect("BostedChristyFitEM-EtaBR", vBReta, false)<7);
722 
723  if (useBRpi1Default || useBRpi2Default || useBRetaDefault)
724  {
725  // use default branching ratios from PDG table
726  // P33(1232)
727  BRpi = 0.;
728  BReta = 0.;
729  BranchingRatios(kPdgP33m1232_DeltaM, brpi, breta);
730  BRpi += brpi;
731  BReta += breta;
732  BranchingRatios(kPdgP33m1232_Delta0, brpi, breta);
733  BRpi += brpi;
734  BReta += breta;
735  BranchingRatios(kPdgP33m1232_DeltaP, brpi, breta);
736  BRpi += brpi;
737  BReta += breta;
738  BranchingRatios(kPdgP33m1232_DeltaPP, brpi, breta);
739  BRpi += brpi;
740  BReta += breta;
741  BRpi /= 4.;
742  BReta /= 4.;
743  fBRp[0][0] = BRpi;
744  fBRp[0][2] = BReta;
745  fBRD[0][0] = BRpi;
746 
747  // S11(1535)
748  BRpi = 0.;
749  BReta = 0.;
750  BranchingRatios(kPdgS11m1535_N0, brpi, breta);
751  BRpi += brpi;
752  BReta += breta;
753  BranchingRatios(kPdgS11m1535_NP, brpi, breta);
754  BRpi += brpi;
755  BReta += breta;
756  BRpi /= 2.;
757  BReta /= 2.;
758  fBRp[1][0] = BRpi;
759  fBRp[1][2] = BReta;
760  fBRD[1][0] = BRpi;
761 
762  // D13(1520)
763  BRpi = 0.;
764  BReta = 0.;
765  BranchingRatios(kPdgD13m1520_N0, brpi, breta);
766  BRpi += brpi;
767  BReta += breta;
768  BranchingRatios(kPdgD13m1520_NP, brpi, breta);
769  BRpi += brpi;
770  BReta += breta;
771  BRpi /= 2.;
772  BReta /= 2.;
773  fBRp[2][0] = BRpi;
774  fBRp[2][2] = BReta;
775  fBRD[2][0] = BRpi;
776 
777  // F15(1680)
778  BRpi = 0.;
779  BReta = 0.;
780  BranchingRatios(kPdgF15m1680_N0, brpi, breta);
781  BRpi += brpi;
782  BReta += breta;
783  BranchingRatios(kPdgF15m1680_NP, brpi, breta);
784  BRpi += brpi;
785  BReta += breta;
786  BRpi /= 2.;
787  BReta /= 2.;
788  fBRp[3][0] = BRpi;
789  fBRp[3][2] = BReta;
790  fBRD[3][0] = BRpi;
791 
792  // S11(1650)
793  BRpi = 0.;
794  BReta = 0.;
795  BranchingRatios(kPdgS11m1650_N0, brpi, breta);
796  BRpi += brpi;
797  BReta += breta;
798  BranchingRatios(kPdgS11m1650_NP, brpi, breta);
799  BRpi += brpi;
800  BReta += breta;
801  BRpi /= 2.;
802  BReta /= 2.;
803  fBRp[4][0] = BRpi;
804  fBRp[4][2] = BReta;
805  fBRD[4][0] = BRpi;
806 
807  // P11(1440)
808  BRpi = 0.;
809  BReta = 0.;
810  BranchingRatios(kPdgP11m1440_N0, brpi, breta);
811  BRpi += brpi;
812  BReta += breta;
813  BranchingRatios(kPdgP11m1440_NP, brpi, breta);
814  BRpi += brpi;
815  BReta += breta;
816  BRpi /= 2.;
817  BReta /= 2.;
818  fBRp[5][0] = BRpi;
819  fBRp[5][2] = BReta;
820  fBRD[5][0] = BRpi;
821 
822  // F37(1950)
823  BRpi = 0.;
824  BReta = 0.;
825  BranchingRatios(kPdgF37m1950_DeltaM , brpi, breta);
826  BRpi += brpi;
827  BReta += breta;
828  BranchingRatios(kPdgF37m1950_Delta0, brpi, breta);
829  BRpi += brpi;
830  BReta += breta;
831  BranchingRatios(kPdgF37m1950_DeltaP, brpi, breta);
832  BRpi += brpi;
833  BReta += breta;
834  BranchingRatios(kPdgF37m1950_DeltaPP, brpi, breta);
835  BRpi += brpi;
836  BReta += breta;
837  BRpi /= 4.;
838  BReta /= 4.;
839  fBRp[6][0] = BRpi;
840  fBRp[6][2] = BReta;
841  fBRD[6][0] = BRpi;
842  }
843  if (!useBRpi1Default)
844  // single pion branching ratios from config file
845  for (int i=0; i<7; i++)
846  fBRp[i][0] = vBRpi1[i];
847  if (!useBRpi2Default)
848  // single pion branching ratios from config file
849  for (int i=0; i<7; i++)
850  fBRD[i][0] = vBRpi2[i];
851  if (!useBRetaDefault)
852  // eta branching ratios from config file
853  for (int i=0; i<7; i++)
854  fBRp[i][2] = vBReta[i];
855 
856  for (int i=0; i<7; i++)
857  fBRD[i][2] = 0.;
858 
859  if (useBRpi1Default || useBRpi2Default)
860  LOG("BostedChristyEMPXSec", pALERT) << "*** Use branching ratios for pion decay from PDG table";
861 
862  if (useBRetaDefault)
863  LOG("BostedChristyEMPXSec", pALERT) << "*** Use branching ratios for eta decay from PDG table";
864 
865  // 2-pion branching ratios
866  for (int i=0;i<7;i++)
867  {
868  fBRp[i][1] = 1.-fBRp[i][0]-fBRp[i][2];
869  fBRD[i][1] = 1.-fBRD[i][0]-fBRD[i][2];
870  }
871 
872  // Meson angular momentum
873  fAngRes[0] = res::OrbitalAngularMom(res::FromPdgCode(kPdgP33m1232_Delta0)); // P33(1232)
874  fAngRes[1] = res::OrbitalAngularMom(res::FromPdgCode(kPdgS11m1535_N0)); // S11(1535)
875  fAngRes[2] = res::OrbitalAngularMom(res::FromPdgCode(kPdgD13m1520_N0)); // D13(1520)
876  fAngRes[3] = res::OrbitalAngularMom(res::FromPdgCode(kPdgF15m1680_N0)); // F15(1680)
877  fAngRes[4] = res::OrbitalAngularMom(res::FromPdgCode(kPdgS11m1650_N0)); // S15(1650)
878  fAngRes[5] = res::OrbitalAngularMom(res::FromPdgCode(kPdgP11m1440_N0)); // P11(1440) roper
879  fAngRes[6] = res::OrbitalAngularMom(res::FromPdgCode(kPdgF37m1950_Delta0)); // F37(1950)
880 
881  std::vector<double> vResMass;
882  // load resonance masses
883  bool useResMassDefault = (GetParamVect("BostedChristyFitEM-ResMass", vResMass, false)<7);
884 
885  if (useResMassDefault)
886  {
887  LOG("BostedChristyEMPXSec", pALERT) << "*** Use resonance masses from PDG table";
888  // Resonance mass
889  fMassRes[0] = res::Mass(res::FromPdgCode(kPdgP33m1232_Delta0)); // P33(1232)
890  fMassRes[1] = res::Mass(res::FromPdgCode(kPdgS11m1535_N0)); // S11(1535)
891  fMassRes[2] = res::Mass(res::FromPdgCode(kPdgD13m1520_N0)); // D13(1520)
892  fMassRes[3] = res::Mass(res::FromPdgCode(kPdgF15m1680_N0)); // F15(1680)
893  fMassRes[4] = res::Mass(res::FromPdgCode(kPdgS11m1650_N0)); // S15(1650)
894  fMassRes[5] = res::Mass(res::FromPdgCode(kPdgP11m1440_N0)); // P11(1440) roper
895  fMassRes[6] = res::Mass(res::FromPdgCode(kPdgF37m1950_Delta0)); // F37(1950)
896  }
897  else
898  // eta branching ratios from config file
899  for (int i=0; i<7; i++)
900  fMassRes[i] = vResMass[i];
901 
902  std::vector<double> vResWidth;
903  // load resonance masses
904  bool useResWidthDefault = (GetParamVect("BostedChristyFitEM-ResWidth", vResWidth, false)<7);
905 
906  if (useResWidthDefault)
907  {
908  LOG("BostedChristyEMPXSec", pALERT) << "*** Use resonance widths from PDG table";
909  // Resonance width
910  fWidthRes[0] = res::Width(res::FromPdgCode(kPdgP33m1232_Delta0)); // P33(1232)
911  fWidthRes[1] = res::Width(res::FromPdgCode(kPdgS11m1535_N0)); // S11(1535)
912  fWidthRes[2] = res::Width(res::FromPdgCode(kPdgD13m1520_N0)); // D13(1520)
913  fWidthRes[3] = res::Width(res::FromPdgCode(kPdgF15m1680_N0)); // F15(1680)
914  fWidthRes[4] = res::Width(res::FromPdgCode(kPdgS11m1650_N0)); // S15(1650)
915  fWidthRes[5] = res::Width(res::FromPdgCode(kPdgP11m1440_N0)); // P11(1440) roper
916  fWidthRes[6] = res::Width(res::FromPdgCode(kPdgF37m1950_Delta0)); // F37(1950)
917  }
918  else
919  // eta branching ratios from config file
920  for (int i=0; i<7; i++)
921  fWidthRes[i] = vResWidth[i];
922 
923  int length;
924 
925  std::vector<double> vRescoef;
926  length = 7;
927  bool isOk = (GetParamVect("BostedChristyFitEM-ResAT0p", vRescoef)>=length);
928  if (!isOk)
929  {
930  LOG("BostedChristyEMPXSec", pFATAL) << "*** Can't find enough proton AT(0)-parameters for xsec^R_T in the config file!";
931  exit(1);
932  }
933  // Ref.1, Table III, AT(0)
934  for (int i=0;i<length;i++)
935  fRescoefTp[i][0] = vRescoef[i];
936 
937  length = 7;
938  isOk = (GetParamVect("BostedChristyFitEM-Resap", vRescoef)>=length);
939  if (!isOk)
940  {
941  LOG("BostedChristyEMPXSec", pFATAL) << "*** Can't find enough proton a-parameters for xsec^R_T in the config file!";
942  exit(1);
943  }
944  // Ref.1, Table III, a
945  for (int i=0;i<length;i++)
946  fRescoefTp[i][1] = vRescoef[i];
947 
948  length = 7;
949  isOk = (GetParamVect("BostedChristyFitEM-Resbp", vRescoef)>=length);
950  if (!isOk)
951  {
952  LOG("BostedChristyEMPXSec", pFATAL) << "*** Can't find enough proton b-parameters parameters for xsec^R_T in the config file!";
953  exit(1);
954  }
955  // Ref.1, Table III, b
956  for (int i=0;i<length;i++)
957  fRescoefTp[i][2] = vRescoef[i];
958 
959  length = 7;
960  isOk = (GetParamVect("BostedChristyFitEM-Rescp", vRescoef)>=length);
961  if (!isOk)
962  {
963  LOG("BostedChristyEMPXSec", pFATAL) << "*** Can't find enough proton c-parameters parameters for xsec^R_T in the config file!";
964  exit(1);
965  }
966  // Ref.1, Table III, c
967  for (int i=0;i<length;i++)
968  fRescoefTp[i][3] = vRescoef[i];
969 
970  length = 7;
971  isOk = (GetParamVect("BostedChristyFitEM-ResAT0D", vRescoef)>=length);
972  if (!isOk)
973  {
974  LOG("BostedChristyEMPXSec", pFATAL) << "*** Can't find enough deuterium AT(0)-parameters for xsec^R_T in the config file!";
975  exit(1);
976  }
977  // Ref.2, Table III, AT(0)
978  for (int i=0;i<length;i++)
979  fRescoefTD[i][0] = vRescoef[i];
980 
981  length = 7;
982  isOk = (GetParamVect("BostedChristyFitEM-ResaD", vRescoef)>=length);
983  if (!isOk)
984  {
985  LOG("BostedChristyEMPXSec", pFATAL) << "*** Can't find enough deuterium a-parameters for xsec^R_T in the config file!";
986  exit(1);
987  }
988  // Ref.2, Table III, a
989  for (int i=0;i<length;i++)
990  fRescoefTD[i][1] = vRescoef[i];
991 
992  length = 7;
993  isOk = (GetParamVect("BostedChristyFitEM-ResbD", vRescoef)>=length);
994  if (!isOk)
995  {
996  LOG("BostedChristyEMPXSec", pFATAL) << "*** Can't find enough deuterium b-parameters parameters for xsec^R_T in the config file!";
997  exit(1);
998  }
999  // Ref.2, Table III, b
1000  for (int i=0;i<length;i++)
1001  fRescoefTD[i][2] = vRescoef[i];
1002 
1003  length = 7;
1004  isOk = (GetParamVect("BostedChristyFitEM-RescD", vRescoef)>=length);
1005  if (!isOk)
1006  {
1007  LOG("BostedChristyEMPXSec", pFATAL) << "*** Can't find enough deuterium c-parameters parameters for xsec^R_T in the config file!";
1008  exit(1);
1009  }
1010  // Ref.2, Table III, c
1011  for (int i=0;i<length;i++)
1012  fRescoefTD[i][3] = vRescoef[i];
1013 
1014  length = 7;
1015  isOk = (GetParamVect("BostedChristyFitEM-ResAL0", vRescoef)>=length);
1016  if (!isOk)
1017  {
1018  LOG("BostedChristyEMPXSec", pFATAL) << "*** Can't find enough proton AL0-parameters parameters for xsec^R_T in the config file!";
1019  exit(1);
1020  }
1021  // Ref.1, Table III, AL(0)
1022  for (int i=0;i<length;i++)
1023  fRescoefL[i][0] = vRescoef[i];
1024 
1025  length = 7;
1026  isOk = (GetParamVect("BostedChristyFitEM-Resd", vRescoef)>=length);
1027  if (!isOk)
1028  {
1029  LOG("BostedChristyEMPXSec", pFATAL) << "*** Can't find enough proton d-parameters parameters for xsec^R_L in the config file!";
1030  exit(1);
1031  }
1032  // Ref.1, Table III, d
1033  for (int i=0;i<length;i++)
1034  fRescoefL[i][1] = vRescoef[i];
1035 
1036  length = 7;
1037  isOk = (GetParamVect("BostedChristyFitEM-Rese", vRescoef)>=length);
1038  if (!isOk)
1039  {
1040  LOG("BostedChristyEMPXSec", pFATAL) << "*** Can't find enough proton e-parameters parameters for xsec^R_L in the config file!";
1041  exit(1);
1042  }
1043  // Ref.1, Table III, e
1044  for (int i=0;i<length;i++)
1045  fRescoefL[i][2] = vRescoef[i];
1046 
1047 
1048  std::vector<double> vNRcoef;
1049  length = 5;
1050  isOk = (GetParamVect("BostedChristyFitEM-NRXSecT1p", vNRcoef)>=length);
1051  if (!isOk)
1052  {
1053  LOG("BostedChristyEMPXSec", pFATAL) << "*** Can't find enough proton bkg parameters for xsec^NR_T in the config file!";
1054  exit(1);
1055  }
1056  // Ref.1, Table IV: \sigma^NR,1_T(0), aT_1, bT_1, cT_1, dT_1
1057  for (int i=0;i<length;i++)
1058  fNRcoefTp[0][i] = vNRcoef[i];
1059 
1060  length = 5;
1061  isOk = (GetParamVect("BostedChristyFitEM-NRXSecT2p", vNRcoef)>=length);
1062  if (!isOk)
1063  {
1064  LOG("BostedChristyEMPXSec", pFATAL) << "*** Can't find enough proton bkg parameters for xsec^NR_T in the config file!";
1065  exit(1);
1066  }
1067  // Ref.1, Table IV: \sigma^NR,2_T(0), aT_2, bT_2, cT_2, dT_2
1068  for (int i=0;i<length;i++)
1069  fNRcoefTp[1][i] = vNRcoef[i];
1070 
1071  length = 5;
1072  isOk = (GetParamVect("BostedChristyFitEM-NRXSecT1D", vNRcoef)>=length);
1073  if (!isOk)
1074  {
1075  LOG("BostedChristyEMPXSec", pFATAL) << "*** Can't find enough deuterium bkg parameters for xsec^NR_T in the config file!";
1076  exit(1);
1077  }
1078  // Ref.2, Table IV: \sigma^NR,1_T(0), aT_1, bT_1, cT_1, dT_1
1079  for (int i=0;i<length;i++)
1080  fNRcoefTD[0][i] = vNRcoef[i];
1081 
1082  length = 5;
1083  isOk = (GetParamVect("BostedChristyFitEM-NRXSecT2D", vNRcoef)>=length);
1084  if (!isOk)
1085  {
1086  LOG("BostedChristyEMPXSec", pFATAL) << "*** Can't find enough deuterium bkg parameters for xsec^NR_T in the config file!";
1087  exit(1);
1088  }
1089  // Ref.2, Table IV: \sigma^NR,2_T(0), aT_2, bT_2, cT_2, dT_2
1090  for (int i=0;i<length;i++)
1091  fNRcoefTD[1][i] = vNRcoef[i];
1092 
1093  length = 6;
1094  isOk = (GetParamVect("BostedChristyFitEM-NRXSecL", vNRcoef)>=length);
1095  if (!isOk)
1096  {
1097  LOG("BostedChristyEMPXSec", pFATAL) << "*** Can't find enough proton bkg parameters for xsec^NR_L in the config file!";
1098  exit(1);
1099  }
1100  // Ref.1, Table IV: \sigma^NR_L, aL, bL, cL, dL, eL
1101  for (int i=0;i<length;i++)
1102  fNRcoefL[i] = vNRcoef[i];
1103 
1104  length = 6;
1105  isOk = (GetParamVect("BostedChristyFitEM-MEC", vNRcoef)>=length);
1106  if (!isOk)
1107  {
1108  LOG("BostedChristyEMPXSec", pFATAL) << "*** Can't find enough parameters for MEC in the config file!";
1109  exit(1);
1110  }
1111  for (int i=0;i<length;i++)
1112  fMECcoef[i] = vNRcoef[i];
1113 
1114  length = 8;
1115  isOk = (GetParamVect("BostedChristyFitEM-MEC2009", vNRcoef)>=length);
1116  if (!isOk)
1117  {
1118  LOG("BostedChristyEMPXSec", pFATAL) << "*** Can't find enough parameters for MEC2009 in the config file!";
1119  exit(1);
1120  }
1121  for (int i=0;i<length;i++)
1122  fMEC2009coef[i] = vNRcoef[i];
1123 
1124  length = 13;
1125  isOk = (GetParamVect("BostedChristyFitEM-Afit", vNRcoef)>=length);
1126  if (!isOk)
1127  {
1128  LOG("BostedChristyEMPXSec", pFATAL) << "*** Can't find enough parameters for nuclei fit (A-fit) in the config file!";
1129  exit(1);
1130  }
1131  for (int i=0;i<length;i++)
1132  fAfitcoef[i] = vNRcoef[i];
1133 
1134 
1135  length = 9;
1136  isOk = (GetParamVect("BostedChristyFitEM-EMCalpha", vNRcoef)>=length);
1137  if (!isOk)
1138  {
1139  LOG("BostedChristyEMPXSec", pFATAL) << "*** Can't find enough alpha coefficients for EMC correction in the config file!";
1140  exit(1);
1141  }
1142  for (int i=0;i<length;i++)
1143  fEMCalpha[i] = vNRcoef[i];
1144 
1145  length = 3;
1146  isOk = (GetParamVect("BostedChristyFitEM-EMCc", vNRcoef)>=length);
1147  if (!isOk)
1148  {
1149  LOG("BostedChristyEMPXSec", pFATAL) << "*** Can't find enough c coefficients for EMC correction in the config file!";
1150  exit(1);
1151  }
1152  for (int i=0;i<length;i++)
1153  fEMCc[i] = vNRcoef[i];
1154 
1155 
1156  std::string keyStart = "BostedChristy-SeparationE@Pdg=";
1157  RgIMap entries = GetConfig().GetItemMap();
1158  for(RgIMap::const_iterator it = entries.begin(); it != entries.end(); ++it)
1159  {
1160  const std::string& key = it->first;
1161  int pdg = 0;
1162  int A = 0;
1163  if (0 == key.compare(0, keyStart.size(), keyStart.c_str()))
1164  {
1165  pdg = atoi(key.c_str() + keyStart.size());
1166  A = pdg::IonPdgCodeToA(pdg);
1167  }
1168  if (0 != pdg && A != 0)
1169  {
1170  std::ostringstream key_ss ;
1171  key_ss << keyStart << pdg;
1172  RgKey rgkey = key_ss.str();
1173  double eb;
1174  GetParam( rgkey, eb) ;
1175  eb = TMath::Max(eb, 0.);
1176  fNucRmvE.insert(map<int,double>::value_type(A,eb));
1177  }
1178  }
1179 
1180  keyStart = "BostedChristy-FermiMomentum@Pdg=";
1181  for(RgIMap::const_iterator it = entries.begin(); it != entries.end(); ++it)
1182  {
1183  const std::string& key = it->first;
1184  int pdg = 0;
1185  int A = 0;
1186  if (0 == key.compare(0, keyStart.size(), keyStart.c_str()))
1187  {
1188  pdg = atoi(key.c_str() + keyStart.size());
1189  A = pdg::IonPdgCodeToA(pdg);
1190  }
1191  if (0 != pdg && A != 0)
1192  {
1193  std::ostringstream key_ss ;
1194  key_ss << keyStart << pdg;
1195  RgKey rgkey = key_ss.str();
1196  double pf;
1197  GetParam( rgkey, pf) ;
1198  pf = TMath::Max(pf, 0.);
1199  fKFTable.insert(map<int,double>::value_type(A,pf));
1200  }
1201  }
1202 
1203  keyStart = "BostedChristy-p18@Pdg=";
1204  for(RgIMap::const_iterator it = entries.begin(); it != entries.end(); ++it)
1205  {
1206  const std::string& key = it->first;
1207  int pdg = 0;
1208  int A = 0;
1209  if (0 == key.compare(0, keyStart.size(), keyStart.c_str()))
1210  {
1211  pdg = atoi(key.c_str() + keyStart.size());
1212  A = pdg::IonPdgCodeToA(pdg);
1213  }
1214  if (0 != pdg && A != 0)
1215  {
1216  std::ostringstream key_ss ;
1217  key_ss << keyStart << pdg;
1218  RgKey rgkey = key_ss.str();
1219  double p18;
1220  GetParam( rgkey, p18) ;
1221  fMEC2009p18.insert(map<int,double>::value_type(A,p18));
1222  }
1223  }
1224 
1225 }
genie::BostedChristyEMPXSec::fMassRes
std::array< double, 7 > fMassRes
resonance mass
Definition: BostedChristyEMPXSec.h:87
genie::ProcessInfo::IsResonant
bool IsResonant(void) const
Definition: ProcessInfo.cxx:99
genie::XSecAlgorithmI
Cross Section Calculation Interface.
Definition: XSecAlgorithmI.h:27
genie::kPdgP33m1232_DeltaPP
const int kPdgP33m1232_DeltaPP
Definition: PDGCodes.h:107
genie::Kinematics::W
double W(bool selected=false) const
Definition: Kinematics.cxx:157
genie::pdg::IsPion
bool IsPion(int pdgc)
Definition: PDGUtils.cxx:326
genie::BostedChristyEMPXSec::Configure
void Configure(const Registry &config)
Definition: BostedChristyEMPXSec.cxx:657
genie::utils::mec::J
double J(double q0, double q3, double Enu, double ml)
Definition: MECUtils.cxx:147
genie::BostedChristyEMPXSec::fRescoefTp
std::array< std::array< double, 4 >, 7 > fRescoefTp
tunable parameters from Ref.1, Table III for resonance
Definition: BostedChristyEMPXSec.h:91
genie::utils::kinematics::Q2
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:1077
genie::Target::HitNucPdg
int HitNucPdg(void) const
Definition: Target.cxx:304
genie::BostedChristyEMPXSec::Integral
double Integral(const Interaction *i) const
Definition: BostedChristyEMPXSec.cxx:613
genie::BostedChristyEMPXSec::fNRcoefTp
std::array< std::array< double, 5 >, 2 > fNRcoefTp
tunable parameters from Ref.1, Table III for nonres bkg
Definition: BostedChristyEMPXSec.h:95
genie::pdg::IsNucleon
bool IsNucleon(int pdgc)
Definition: PDGUtils.cxx:346
genie::Target::A
int A(void) const
Definition: Target.h:70
genie::BostedChristyEMPXSec::fEMCalpha
std::array< double, 9 > fEMCalpha
tunable parameters for EMC fit
Definition: BostedChristyEMPXSec.h:102
genie::kPdgF15m1680_NP
const int kPdgF15m1680_NP
Definition: PDGCodes.h:135
genie::BostedChristyEMPXSec::sigmaR
double sigmaR(int, double, double, bool) const
Definition: BostedChristyEMPXSec.cxx:53
pFATAL
#define pFATAL
Definition: Messenger.h:56
genie::BostedChristyEMPXSec::fMEC2009coef
std::array< double, 8 > fMEC2009coef
tunable parameters for MEC2009 function
Definition: BostedChristyEMPXSec.h:99
genie::Algorithm::GetParamVect
int GetParamVect(const std::string &comm_name, std::vector< T > &v, bool is_top_call=true) const
Handle to load vectors of parameters.
genie::pdg::IonPdgCodeToA
int IonPdgCodeToA(int pdgc)
Definition: PDGUtils.cxx:63
genie::Kinematics
Generated/set kinematical variables for an event.
Definition: Kinematics.h:39
genie::pdg::IsChargedLepton
bool IsChargedLepton(int pdgc)
Definition: PDGUtils.cxx:101
genie::BostedChristyEMPXSec::fEMCc
std::array< double, 3 > fEMCc
tunable parameters for EMC fit
Definition: BostedChristyEMPXSec.h:103
genie::utils::res::Mass
double Mass(Resonance_t res)
resonance mass (GeV)
Definition: BaryonResUtils.cxx:436
genie::kPdgS11m1650_N0
const int kPdgS11m1650_N0
Definition: PDGCodes.h:112
genie::BostedChristyEMPXSec::fBRD
std::array< std::array< double, 3 >, 7 > fBRD
branching ratios of resonances for deterium fit
Definition: BostedChristyEMPXSec.h:83
genie::utils::res::Width
double Width(Resonance_t res)
resonance width (GeV)
Definition: BaryonResUtils.cxx:463
genie::kPdgF37m1950_DeltaM
const int kPdgF37m1950_DeltaM
Definition: PDGCodes.h:148
genie::kPdgF37m1950_DeltaP
const int kPdgF37m1950_DeltaP
Definition: PDGCodes.h:150
genie::KinePhaseSpace_t
enum genie::EKinePhaseSpace KinePhaseSpace_t
genie::units::ub
static constexpr double ub
Definition: Units.h:80
genie::kPdgS11m1535_N0
const int kPdgS11m1535_N0
Definition: PDGCodes.h:108
genie::Algorithm::GetConfig
virtual const Registry & GetConfig(void) const
Definition: Algorithm.cxx:246
genie::BostedChristyEMPXSec::fMECcoef
std::array< double, 6 > fMECcoef
tunable parameters for Eqs.(20), (21) Ref.2
Definition: BostedChristyEMPXSec.h:98
genie::kPdgP33m1232_DeltaP
const int kPdgP33m1232_DeltaP
Definition: PDGCodes.h:106
genie::kPdgP33m1232_DeltaM
const int kPdgP33m1232_DeltaM
Definition: PDGCodes.h:104
genie::BostedChristyEMPXSec::FermiSmearingA
void FermiSmearingA(double, double, double, double, double &, double &, double &, double &) const
Definition: BostedChristyEMPXSec.cxx:199
genie::utils::kinematics::W
double W(const Interaction *const i)
Definition: KineUtils.cxx:1101
genie::pdg::IsNeutron
bool IsNeutron(int pdgc)
Definition: PDGUtils.cxx:341
pALERT
#define pALERT
Definition: Messenger.h:57
genie::Interaction
Summary information for an interaction.
Definition: Interaction.h:56
genie::BostedChristyEMPXSec::fAfitcoef
std::array< double, 13 > fAfitcoef
tunable parameters for nuclei fit
Definition: BostedChristyEMPXSec.h:100
genie::kPdgF37m1950_Delta0
const int kPdgF37m1950_Delta0
Definition: PDGCodes.h:149
genie::pdg::IsProton
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:336
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
genie::Registry::GetItemMap
const RgIMap & GetItemMap(void) const
Definition: Registry.h:161
genie::units::A
static constexpr double A
Definition: Units.h:74
genie::ProcessInfo
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition: ProcessInfo.h:46
genie::Target
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:40
genie::Algorithm::Configure
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
genie::InitialState::ProbePdg
int ProbePdg(void) const
Definition: InitialState.h:64
genie::BostedChristyEMPXSec::BranchingRatios
void BranchingRatios(int, double &, double &) const
Definition: BostedChristyEMPXSec.cxx:669
genie::kPdgEta
const int kPdgEta
Definition: PDGCodes.h:161
genie::kPdgPi0
const int kPdgPi0
Definition: PDGCodes.h:160
genie::utils::res::OrbitalAngularMom
int OrbitalAngularMom(Resonance_t res)
orbital angular momentum
Definition: BaryonResUtils.cxx:531
genie::Target::Z
int Z(void) const
Definition: Target.h:68
genie::kPdgF37m1950_DeltaPP
const int kPdgF37m1950_DeltaPP
Definition: PDGCodes.h:151
genie::BostedChristyEMPXSec::fBRp
std::array< std::array< double, 3 >, 7 > fBRp
branching ratios of resonances for proton fit
Definition: BostedChristyEMPXSec.h:82
genie::kPdgD13m1520_NP
const int kPdgD13m1520_NP
Definition: PDGCodes.h:111
genie::kPdgF15m1680_N0
const int kPdgF15m1680_N0
Definition: PDGCodes.h:134
genie::BostedChristyEMPXSec::fWidthRes
std::array< double, 7 > fWidthRes
resonance width
Definition: BostedChristyEMPXSec.h:89
genie::kPdgS11m1650_NP
const int kPdgS11m1650_NP
Definition: PDGCodes.h:113
genie::utils::res::FromPdgCode
Resonance_t FromPdgCode(int pdgc)
PDG code -&gt; resonance id.
Definition: BaryonResUtils.cxx:84
genie::kPdgP33m1232_Delta0
const int kPdgP33m1232_Delta0
Definition: PDGCodes.h:105
genie::kPdgP11m1440_N0
const int kPdgP11m1440_N0
Definition: PDGCodes.h:126
genie::BostedChristyEMPXSec::fNRcoefTD
std::array< std::array< double, 5 >, 2 > fNRcoefTD
tunable parameters from Ref.1, Table IV for nonres bkg
Definition: BostedChristyEMPXSec.h:96
genie::ProcessInfo::IsEM
bool IsEM(void) const
Definition: ProcessInfo.cxx:193
genie::kPdgD13m1520_N0
const int kPdgD13m1520_N0
Definition: PDGCodes.h:110
genie::BostedChristyEMPXSec::fXSecIntegrator
const XSecIntegratorI * fXSecIntegrator
Definition: BostedChristyEMPXSec.h:109
genie::kPdgP11m1440_NP
const int kPdgP11m1440_NP
Definition: PDGCodes.h:127
genie::PDGLibrary::Instance
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:68
genie::BostedChristyEMPXSec::fUseMEC
bool fUseMEC
account for MEC contribution?
Definition: BostedChristyEMPXSec.h:70
genie::PDGLibrary
Singleton class to load &amp; serve a TDatabasePDG.
Definition: PDGLibrary.h:35
RgKey
string RgKey
Definition: RegistryItemTypeDef.h:33
genie::Registry
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:65
genie::BostedChristyEMPXSec::fRescoefTD
std::array< std::array< double, 4 >, 7 > fRescoefTD
tunable parameters from Ref.2, Table III for resonance
Definition: BostedChristyEMPXSec.h:92
genie::utils::kinematics::Jacobian
double Jacobian(const Interaction *const i, KinePhaseSpace_t f, KinePhaseSpace_t t)
Definition: KineUtils.cxx:130
genie::BostedChristyEMPXSec::FitEMC
double FitEMC(double, int) const
Definition: BostedChristyEMPXSec.cxx:582
genie::BostedChristyEMPXSec::fNRcoefL
std::array< double, 6 > fNRcoefL
tunable parameters from Ref.1, Table III for nonres bkg
Definition: BostedChristyEMPXSec.h:97
genie::Interaction::InitState
const InitialState & InitState(void) const
Definition: Interaction.h:69
genie::Interaction::ProcInfo
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
genie::kPdgProton
const int kPdgProton
Definition: PDGCodes.h:81
genie::PDGLibrary::Find
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:86
genie::BostedChristyEMPXSec::fAngRes
std::array< int, 7 > fAngRes
resonance angular momentum
Definition: BostedChristyEMPXSec.h:85
genie::BostedChristyEMPXSec::MEC2009
double MEC2009(int, double, double) const
Definition: BostedChristyEMPXSec.cxx:407
genie::Kinematics::Q2
double Q2(bool selected=false) const
Definition: Kinematics.cxx:125
genie::Algorithm::GetParamDef
bool GetParamDef(const RgKey &name, T &p, const T &def) const
genie::Algorithm::GetParam
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
genie::InitialState::Tgt
const Target & Tgt(void) const
Definition: InitialState.h:66
genie::BostedChristyEMPXSec::FermiSmearingD
void FermiSmearingD(double, double, double &, double &, double &, double &, bool) const
Definition: BostedChristyEMPXSec.cxx:269
genie::kPdgS11m1535_NP
const int kPdgS11m1535_NP
Definition: PDGCodes.h:109
genie::BostedChristyEMPXSec::ValidKinematics
bool ValidKinematics(const Interaction *i) const
Is the input kinematical point a physically allowed one?
Definition: BostedChristyEMPXSec.cxx:644
genie::kPdgTgtDeuterium
const int kPdgTgtDeuterium
Definition: PDGCodes.h:201
genie::XSecIntegratorI::Integrate
virtual double Integrate(const XSecAlgorithmI *model, const Interaction *interaction) const =0
genie::InitialState::ProbeE
double ProbeE(RefFrame_t rf) const
Definition: InitialState.cxx:384
PDGCodes.h
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
genie::BostedChristyEMPXSec::fRescoefL
std::array< std::array< double, 3 >, 7 > fRescoefL
tunable parameters from Ref.1, Table III for resonance
Definition: BostedChristyEMPXSec.h:93
genie::BostedChristyEMPXSec::ValidProcess
bool ValidProcess(const Interaction *i) const
Can this cross section algorithm handle the input process?
Definition: BostedChristyEMPXSec.cxx:619
genie::kISkipProcessChk
const UInt_t kISkipProcessChk
if set, skip process validity checks
Definition: Interaction.h:47
genie::BostedChristyEMPXSec::sigmaNR
double sigmaNR(int, double, double, bool) const
Definition: BostedChristyEMPXSec.cxx:155
genie::BostedChristyEMPXSec::XSec
double XSec(const Interaction *i, KinePhaseSpace_t k) const
Compute the cross section for the input interaction.
Definition: BostedChristyEMPXSec.cxx:439
genie::InitialState
Initial State information.
Definition: InitialState.h:48
genie::RgIMap
map< RgKey, RegistryItemI * > RgIMap
Definition: Registry.h:45
Generated on Wed Jul 17 2024 14:46:38 for GENIEGenerator by   doxygen 1.8.5