GENIEGenerator
Loading...
Searching...
No Matches
BostedChristyEMPXSec.cxx
Go to the documentation of this file.
1//_________________________________________________________________________
2/*
3 Copyright (c) 2003-2025, 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
30using namespace genie;
31using namespace genie::constants;
32using namespace genie::utils;
33
34//____________________________________________________________________________
35BostedChristyEMPXSec::BostedChristyEMPXSec() :
36XSecAlgorithmI("genie::BostedChristyEMPXSec")
37{
38}
39//____________________________________________________________________________
40BostedChristyEMPXSec::BostedChristyEMPXSec(string config) :
41XSecAlgorithmI("genie::BostedChristyEMPXSec", config)
42{
43}
44//____________________________________________________________________________
49//____________________________________________________________________________
50// resonance cross section fit of electron-proton and electron-deuterium scattering data
51// sf=0 \sigma_T
52// sf=1 \sigma_L
53double 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
155double 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
199void 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
269void 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//____________________________________________________________________________
407double 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//____________________________________________________________________________
439double 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//____________________________________________________________________________
582double 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//____________________________________________________________________________
613double BostedChristyEMPXSec::Integral(const Interaction * interaction) const
614{
615 double xsec = fXSecIntegrator->Integrate(this,interaction);
616 return xsec;
617}
618//____________________________________________________________________________
619bool 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//____________________________________________________________________________
644bool 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//____________________________________________________________________________
662//____________________________________________________________________________
668//____________________________________________________________________________
669void 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//____________________________________________________________________________
694void 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}
pALERT
#define pALERT
Definition Messenger.h:57
pFATAL
#define pFATAL
Definition Messenger.h:56
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE,...
Definition Messenger.h:96
PDGCodes.h
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils.
RgKey
string RgKey
Definition RegistryItemTypeDef.h:33
genie::Algorithm::GetConfig
virtual const Registry & GetConfig(void) const
Definition Algorithm.cxx:246
genie::Algorithm::GetParam
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
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::Algorithm::Configure
virtual void Configure(const Registry &config)
Definition Algorithm.cxx:62
genie::Algorithm::GetParamDef
bool GetParamDef(const RgKey &name, T &p, const T &def) const
genie::BostedChristyEMPXSec::FitEMC
double FitEMC(double, int) const
Definition BostedChristyEMPXSec.cxx:582
genie::BostedChristyEMPXSec::fRescoefTp
std::array< std::array< double, 4 >, 7 > fRescoefTp
tunable parameters from Ref.1, Table III for resonance \sigma_T
Definition BostedChristyEMPXSec.h:91
genie::BostedChristyEMPXSec::ValidProcess
bool ValidProcess(const Interaction *i) const
Can this cross section algorithm handle the input process?
Definition BostedChristyEMPXSec.cxx:619
genie::BostedChristyEMPXSec::fNRcoefTp
std::array< std::array< double, 5 >, 2 > fNRcoefTp
tunable parameters from Ref.1, Table III for nonres bkg \sigma_T
Definition BostedChristyEMPXSec.h:95
genie::BostedChristyEMPXSec::MEC2009
double MEC2009(int, double, double) const
Definition BostedChristyEMPXSec.cxx:407
genie::BostedChristyEMPXSec::Integral
double Integral(const Interaction *i) const
Definition BostedChristyEMPXSec.cxx:613
genie::BostedChristyEMPXSec::fXSecIntegrator
const XSecIntegratorI * fXSecIntegrator
Definition BostedChristyEMPXSec.h:109
genie::BostedChristyEMPXSec::fMassRes
std::array< double, 7 > fMassRes
resonance mass
Definition BostedChristyEMPXSec.h:87
genie::BostedChristyEMPXSec::fEMCalpha
std::array< double, 9 > fEMCalpha
tunable parameters for EMC fit
Definition BostedChristyEMPXSec.h:102
genie::BostedChristyEMPXSec::fAfitcoef
std::array< double, 13 > fAfitcoef
tunable parameters for nuclei fit
Definition BostedChristyEMPXSec.h:100
genie::BostedChristyEMPXSec::fNRcoefTD
std::array< std::array< double, 5 >, 2 > fNRcoefTD
tunable parameters from Ref.1, Table IV for nonres bkg \sigma_T
Definition BostedChristyEMPXSec.h:96
genie::BostedChristyEMPXSec::fWidthRes
std::array< double, 7 > fWidthRes
resonance width
Definition BostedChristyEMPXSec.h:89
genie::BostedChristyEMPXSec::fRescoefTD
std::array< std::array< double, 4 >, 7 > fRescoefTD
tunable parameters from Ref.2, Table III for resonance \sigma_T
Definition BostedChristyEMPXSec.h:92
genie::BostedChristyEMPXSec::fRescoefL
std::array< std::array< double, 3 >, 7 > fRescoefL
tunable parameters from Ref.1, Table III for resonance \sigma_L
Definition BostedChristyEMPXSec.h:93
genie::BostedChristyEMPXSec::BranchingRatios
void BranchingRatios(int, double &, double &) const
Definition BostedChristyEMPXSec.cxx:669
genie::BostedChristyEMPXSec::Configure
void Configure(const Registry &config)
Definition BostedChristyEMPXSec.cxx:657
genie::BostedChristyEMPXSec::FermiSmearingD
void FermiSmearingD(double, double, double &, double &, double &, double &, bool) const
Definition BostedChristyEMPXSec.cxx:269
genie::BostedChristyEMPXSec::fUseMEC
bool fUseMEC
account for MEC contribution?
Definition BostedChristyEMPXSec.h:70
genie::BostedChristyEMPXSec::fEMCc
std::array< double, 3 > fEMCc
tunable parameters for EMC fit
Definition BostedChristyEMPXSec.h:103
genie::BostedChristyEMPXSec::sigmaNR
double sigmaNR(int, double, double, bool) const
Definition BostedChristyEMPXSec.cxx:155
genie::BostedChristyEMPXSec::fNRcoefL
std::array< double, 6 > fNRcoefL
tunable parameters from Ref.1, Table III for nonres bkg \sigma_L
Definition BostedChristyEMPXSec.h:97
genie::BostedChristyEMPXSec::fBRp
std::array< std::array< double, 3 >, 7 > fBRp
branching ratios of resonances for proton fit
Definition BostedChristyEMPXSec.h:82
genie::BostedChristyEMPXSec::fMECcoef
std::array< double, 6 > fMECcoef
tunable parameters for Eqs.(20), (21) Ref.2
Definition BostedChristyEMPXSec.h:98
genie::BostedChristyEMPXSec::fBRD
std::array< std::array< double, 3 >, 7 > fBRD
branching ratios of resonances for deterium fit
Definition BostedChristyEMPXSec.h:83
genie::BostedChristyEMPXSec::FermiSmearingA
void FermiSmearingA(double, double, double, double, double &, double &, double &, double &) const
Definition BostedChristyEMPXSec.cxx:199
genie::BostedChristyEMPXSec::fAngRes
std::array< int, 7 > fAngRes
resonance angular momentum
Definition BostedChristyEMPXSec.h:85
genie::BostedChristyEMPXSec::fMEC2009coef
std::array< double, 8 > fMEC2009coef
tunable parameters for MEC2009 function
Definition BostedChristyEMPXSec.h:99
genie::BostedChristyEMPXSec::ValidKinematics
bool ValidKinematics(const Interaction *i) const
Is the input kinematical point a physically allowed one?
Definition BostedChristyEMPXSec.cxx:644
genie::BostedChristyEMPXSec::sigmaR
double sigmaR(int, double, double, bool) const
Definition BostedChristyEMPXSec.cxx:53
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::InitialState::Tgt
const Target & Tgt(void) const
Definition InitialState.h:66
genie::InitialState::ProbePdg
int ProbePdg(void) const
Definition InitialState.h:64
genie::InitialState::ProbeE
double ProbeE(RefFrame_t rf) const
Definition InitialState.cxx:384
genie::Interaction
Summary information for an interaction.
Definition Interaction.h:56
genie::Interaction::ProcInfo
const ProcessInfo & ProcInfo(void) const
Definition Interaction.h:70
genie::Interaction::InitState
const InitialState & InitState(void) const
Definition Interaction.h:69
genie::Kinematics
Generated/set kinematical variables for an event.
Definition Kinematics.h:39
genie::PDGLibrary
Singleton class to load & serve a TDatabasePDG.
Definition PDGLibrary.h:36
genie::PDGLibrary::Instance
static PDGLibrary * Instance(void)
Definition PDGLibrary.cxx:68
genie::PDGLibrary::Find
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition PDGLibrary.cxx:86
genie::ProcessInfo
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition ProcessInfo.h:46
genie::ProcessInfo::IsEM
bool IsEM(void) const
Definition ProcessInfo.cxx:193
genie::ProcessInfo::IsResonant
bool IsResonant(void) const
Definition ProcessInfo.cxx:99
genie::Registry
A registry. Provides the container for algorithm configuration parameters.
Definition Registry.h:65
genie::Registry::GetItemMap
const RgIMap & GetItemMap(void) const
Definition Registry.h:161
genie::Target
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition Target.h:40
genie::Target::HitNucPdg
int HitNucPdg(void) const
Definition Target.cxx:304
genie::Target::Z
int Z(void) const
Definition Target.h:68
genie::Target::A
int A(void) const
Definition Target.h:70
genie::constants
Basic constants.
genie::pdg
Utilities for improving the code readability when using PDG codes.
genie::pdg::IsProton
bool IsProton(int pdgc)
Definition PDGUtils.cxx:336
genie::pdg::IsNucleon
bool IsNucleon(int pdgc)
Definition PDGUtils.cxx:346
genie::pdg::IsChargedLepton
bool IsChargedLepton(int pdgc)
Definition PDGUtils.cxx:101
genie::pdg::IsNeutron
bool IsNeutron(int pdgc)
Definition PDGUtils.cxx:341
genie::pdg::IsPion
bool IsPion(int pdgc)
Definition PDGUtils.cxx:326
genie::pdg::IonPdgCodeToA
int IonPdgCodeToA(int pdgc)
Definition PDGUtils.cxx:63
genie::units::ub
static constexpr double ub
Definition Units.h:80
genie::utils::config
Simple functions for loading and reading nucleus dependent keys from config files.
Definition ConfigIsotopeMapUtils.cxx:21
genie::utils::kinematics
Kinematical utilities.
genie::utils::kinematics::Jacobian
double Jacobian(const Interaction *const i, KinePhaseSpace_t f, KinePhaseSpace_t t)
Definition KineUtils.cxx:130
genie::utils::res::FromPdgCode
Resonance_t FromPdgCode(int pdgc)
PDG code -> resonance id.
Definition BaryonResUtils.cxx:84
genie::utils::res::Width
double Width(Resonance_t res)
resonance width (GeV)
Definition BaryonResUtils.cxx:463
genie::utils::res::Mass
double Mass(Resonance_t res)
resonance mass (GeV)
Definition BaryonResUtils.cxx:436
genie::utils::res::OrbitalAngularMom
int OrbitalAngularMom(Resonance_t res)
orbital angular momentum
Definition BaryonResUtils.cxx:531
genie::utils
Root of GENIE utility namespaces.
genie
THE MAIN GENIE PROJECT NAMESPACE
Definition AlgCmp.h:25
genie::kPdgEta
const int kPdgEta
Definition PDGCodes.h:161
genie::kPdgD13m1520_NP
const int kPdgD13m1520_NP
Definition PDGCodes.h:111
genie::kPdgF37m1950_Delta0
const int kPdgF37m1950_Delta0
Definition PDGCodes.h:149
genie::kPdgD13m1520_N0
const int kPdgD13m1520_N0
Definition PDGCodes.h:110
genie::kPdgP33m1232_Delta0
const int kPdgP33m1232_Delta0
Definition PDGCodes.h:105
genie::RgIMap
map< RgKey, RegistryItemI * > RgIMap
Definition Registry.h:45
genie::kPdgProton
const int kPdgProton
Definition PDGCodes.h:81
genie::kPdgPi0
const int kPdgPi0
Definition PDGCodes.h:160
genie::kPdgTgtDeuterium
const int kPdgTgtDeuterium
Definition PDGCodes.h:201
genie::kPdgF15m1680_N0
const int kPdgF15m1680_N0
Definition PDGCodes.h:134
genie::kPdgP33m1232_DeltaPP
const int kPdgP33m1232_DeltaPP
Definition PDGCodes.h:107
genie::kPdgP11m1440_N0
const int kPdgP11m1440_N0
Definition PDGCodes.h:126
genie::kPdgS11m1535_NP
const int kPdgS11m1535_NP
Definition PDGCodes.h:109
genie::kPdgF37m1950_DeltaM
const int kPdgF37m1950_DeltaM
Definition PDGCodes.h:148
genie::kPdgS11m1650_N0
const int kPdgS11m1650_N0
Definition PDGCodes.h:112
genie::kPdgF37m1950_DeltaP
const int kPdgF37m1950_DeltaP
Definition PDGCodes.h:150
genie::kPdgF37m1950_DeltaPP
const int kPdgF37m1950_DeltaPP
Definition PDGCodes.h:151
genie::KinePhaseSpace_t
enum genie::EKinePhaseSpace KinePhaseSpace_t
genie::kPdgS11m1535_N0
const int kPdgS11m1535_N0
Definition PDGCodes.h:108
genie::kPdgS11m1650_NP
const int kPdgS11m1650_NP
Definition PDGCodes.h:113
genie::kPdgP11m1440_NP
const int kPdgP11m1440_NP
Definition PDGCodes.h:127
genie::kPdgP33m1232_DeltaM
const int kPdgP33m1232_DeltaM
Definition PDGCodes.h:104
genie::kRfHitNucRest
@ kRfHitNucRest
Definition RefFrame.h:30
genie::kISkipProcessChk
const UInt_t kISkipProcessChk
if set, skip process validity checks
Definition Interaction.h:47
genie::kPdgP33m1232_DeltaP
const int kPdgP33m1232_DeltaP
Definition PDGCodes.h:106
genie::kPdgF15m1680_NP
const int kPdgF15m1680_NP
Definition PDGCodes.h:135
Generated on for GENIEGenerator by doxygen 1.14.0