GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
TabulatedLabFrameHadronTensor.cxx
Go to the documentation of this file.
1 // standard library includes
2 #include <cmath>
3 #include <fstream>
4 
5 // GENIE includes
11 
12 // For retrieval of CKM-Vud
15 
16 namespace {
17  /// Enumerated type that represents the format used to read in
18  /// grid points from the hadron tensor data file
20  /// A starting value and a step size are used to define a regular grid
21  kStartAndStep = 0,
22  /// An explicit table of grid points is used to define an irregular grid
23  kExplicitValues = 1,
24  /// Dummy value used for error checking
25  kHadronTensorGridFlag_COUNT = 2
26  };
27 
28  /// Definition of sqrt() that returns zero if the argument is negative.
29  /// Used to prevent spurious NaNs due to numerical roundoff.
30  double real_sqrt(double x) {
31  if (x < 0.) return 0.;
32  else return std::sqrt(x);
33  }
34 }
35 
37  const std::string& table_file_name)
38  : fGrid(&fq0Points, &fqmagPoints, &fEntries)
39 {
40  // Read in the table
41  std::ifstream in_file( table_file_name.c_str() );
42 
43 
44  // Skip the initial comment line
45  std::string dummy;
46  std::getline(in_file, dummy);
47 
48  /// \todo Add error checks
49  std::string type_name;
50  int Z, A, num_q0, num_q_mag;
51 
52  /// \todo Use type name
53  in_file >> Z >> A >> type_name >> num_q0 >> num_q_mag;
54 
55  int q0_flag;
56  in_file >> q0_flag;
57  read1DGridValues(num_q0, q0_flag, in_file, fq0Points);
58 
59  int q_mag_flag;
60  in_file >> q_mag_flag;
61  read1DGridValues(num_q_mag, q_mag_flag, in_file, fqmagPoints);
62 
64 
65  std::cout<< "TabulatedLabFrameHadronTensor: reading hadron tensor table" << std::endl;
66  std::cout<< "table_file_name: " << table_file_name << std::endl;
67  std::cout<< "Z: " << Z << std::endl;
68  std::cout<< "A: " << A << std::endl;
69  std::cout<< "num_q0: " << num_q0 << std::endl;
70  std::cout<< "num_q_mag: " << num_q_mag << std::endl;
71 
72  double W00, ReW0z, Wxx, ImWxy, Wzz;
73  int lineCount=1;
74 
75  for (long j = 0; j < num_q0; ++j) {
76  for (long k = 0; k < num_q_mag; ++k) {
77 
78  fEntries.push_back( TableEntry() );
79  TableEntry& entry = fEntries.back();
80 
81  // in_file >> entry.W00 >> entry.ReW0z >> entry.Wxx
82  // >> entry.ImWxy >> entry.Wzz;
83 
84  in_file >> W00 >> ReW0z >> Wxx >> ImWxy >> Wzz;
85 
86  entry.W00 = W00;
87  entry.ReW0z= ReW0z;
88  entry.Wxx = Wxx;
89  entry.ImWxy= ImWxy;
90  entry.Wzz = Wzz;
91 
92  //if(W00 > 0.0) {
93  // std::cout<< "line count is" << lineCount << std::endl;
94  // std::cout<< "Entry: " << entry.W00 << " " << entry.ReW0z << " " << entry.Wxx << " " << entry.ImWxy << " " << entry.Wzz << std::endl;
95  // std::cout<< "File: " << W00 << " " << ReW0z << " " << Wxx << " " << ImWxy << " " << Wzz << std::endl;
96  //}
97  lineCount++;
98  }
99  }
100 }
101 
103 {
104 }
105 
107  double q0, double q_mag) const
108 {
109  TableEntry entry = fGrid.interpolate(q0, q_mag);
110  return std::complex<double>(entry.W00, 0.);
111 }
112 
114  double q0, double q_mag) const
115 {
116  TableEntry entry = fGrid.interpolate(q0, q_mag);
117  // Currently only the real part of W0z is tabulated
118  /// \todo Think about adding the imaginary part even though it's not needed
119  /// for the cross section calculation
120  return std::complex<double>(entry.ReW0z, 0.);
121 }
122 
124  double q0, double q_mag) const
125 {
126  TableEntry entry = fGrid.interpolate(q0, q_mag);
127  return std::complex<double>(entry.Wxx, 0.);
128 }
129 
131  double q0, double q_mag) const
132 {
133  TableEntry entry = fGrid.interpolate(q0, q_mag);
134  // The Wxy element is purely imaginary
135  return std::complex<double>(0., entry.ImWxy);
136 }
137 
139  double q0, double q_mag) const
140 {
141  TableEntry entry = fGrid.interpolate(q0, q_mag);
142  // The Wxy element is purely imaginary
143  return std::complex<double>(entry.Wzz, 0.);
144 }
145 
147  double q_mag, double Mi) const
148 {
149  TableEntry entry = fGrid.interpolate(q0, q_mag);
150  return W1(q0, q_mag, entry) / Mi;
151 }
152 
154  double q_mag, double Mi) const
155 {
156  TableEntry entry = fGrid.interpolate(q0, q_mag);
157  return W2(q0, q_mag, entry) / Mi;
158 }
159 
161  double q_mag, double /*Mi*/) const
162 {
163  TableEntry entry = fGrid.interpolate(q0, q_mag);
164  return W3(q0, q_mag, entry);
165 }
166 
168  double q_mag, double Mi) const
169 {
170  TableEntry entry = fGrid.interpolate(q0, q_mag);
171  return W4(q0, q_mag, entry) * Mi;
172 }
173 
175  double q_mag, double /*Mi*/) const
176 {
177  TableEntry entry = fGrid.interpolate(q0, q_mag);
178  return W5(q0, q_mag, entry);
179 }
180 
182  double /* q_mag */, double /* Mi */) const
183 {
184  return 0.;
185 }
186 
188  int flag, std::ifstream& in_file, std::vector<double>& vec_to_fill)
189 {
190  vec_to_fill.clear();
191 
192  if (flag >= kHadronTensorGridFlag_COUNT) {
193  LOG("TabulatedLabFrameHadronTensor", pERROR)
194  << "Invalid hadron tensor grid flag value \"" << flag
195  << "\" encountered.";
196  return;
197  }
198  else if (flag == kStartAndStep) {
199  double start_val, step;
200  in_file >> start_val >> step;
201  for (int k = 0; k < num_points; ++k)
202  vec_to_fill.push_back( start_val + (k * step) );
203  }
204  else if (flag == kExplicitValues) {
205  double val;
206  for (int k = 0; k < num_points; ++k) {
207  in_file >> val;
208  vec_to_fill.push_back( val );
209  }
210  }
211 
212 }
213 
215  double /*q_mag*/, const TableEntry& entry) const
216 {
217  return entry.Wxx / 2.;
218 }
219 
220 double genie::TabulatedLabFrameHadronTensor::W2(double q0, double q_mag,
221  const TableEntry& entry) const
222 {
223  double temp_1 = ( std::pow(q0, 2) / std::pow(q_mag, 2) )
224  * (entry.Wzz - entry.Wxx);
225  double temp_2 = 2. * q0 * entry.ReW0z / q_mag;
226  return ( entry.W00 + entry.Wxx + temp_1 - temp_2 ) / 2.;
227 }
228 
229 double genie::TabulatedLabFrameHadronTensor::W3(double /*q0*/, double q_mag,
230  const TableEntry& entry) const
231 {
232  return ( entry.ImWxy / q_mag );
233 }
234 
235 double genie::TabulatedLabFrameHadronTensor::W4(double /*q0*/, double q_mag,
236  const TableEntry& entry) const
237 {
238  return ( entry.Wzz - entry.Wxx ) / ( 2. * std::pow(q_mag, 2) );
239 }
240 
241 double genie::TabulatedLabFrameHadronTensor::W5(double q0, double q_mag,
242  const TableEntry& entry) const
243 {
244  return ( entry.ReW0z - q0 * (entry.Wzz - entry.Wxx) / q_mag ) / q_mag;
245 }
246 
248  double /*q_mag*/, const TableEntry& /*entry*/) const
249 {
250  // Currently, \Im W^{0z} has not been tabulated, so we can't really
251  // evaluate W6. This isn't a huge problem, however, because W6 doesn't
252  // contribute to neutrino cross sections.
253  return 0.;
254 }
255 
257  const Interaction* interaction, double Q_value) const
258 {
259  // Don't do anything if you've been handed a nullptr
260  if ( !interaction ) return 0.;
261 
262  int probe_pdg = interaction->InitState().ProbePdg();
263  double E_probe = interaction->InitState().ProbeE(kRfLab);
264  double m_probe = interaction->InitState().Probe()->Mass();
265  double Tl = interaction->Kine().GetKV(kKVTl);
266  double cos_l = interaction->Kine().GetKV(kKVctl);
267  double ml = interaction->FSPrimLepton()->Mass();
268 
269  return dSigma_dT_dCosTheta(probe_pdg, E_probe, m_probe, Tl, cos_l, ml,
270  Q_value);
271 }
272 
274  double E_probe, double m_probe, double Tl, double cos_l, double ml,
275  double Q_value) const
276 {
277  // dSigma_dT_dCosTheta in GeV^(-3)
278  double xsec = 0.;
279 
280  /// \todo Add check if in grid. If not, return 0
281 
282  // Final state lepton total energy
283  double El = Tl + ml;
284 
285  // Energy transfer (uncorrected)
286  double q0 = E_probe - El;
287 
288  // The corrected energy transfer takes into account the binding
289  // energy of the struck nucleon(s)
290  double q0_corrected = q0 - Q_value;
291 
292  // Magnitude of the initial state lepton 3-momentum
293  double k_initial = real_sqrt( std::pow(E_probe, 2) - std::pow(m_probe, 2) );
294 
295  // Magnitude of the final state lepton 3-momentum
296  double k_final = real_sqrt( std::pow(Tl, 2) + 2*ml*Tl );
297 
298  // Magnitude of the 3-momentum transfer
299  double q_mag2 = std::pow(k_initial, 2) + std::pow(k_final, 2)
300  - 2.*k_initial*k_final*cos_l;
301  double q_mag = real_sqrt( q_mag2 );
302 
303  // Find the appropriate values of the hadron tensor elements for the
304  // given combination of q0_corrected and q_mag
305  TableEntry entry = fGrid.interpolate(q0_corrected, q_mag);
306 
307  // The half-angle formulas come in handy here. See, e.g.,
308  // http://mathworld.wolfram.com/Half-AngleFormulas.html
309  double s2_half = (1. - cos_l) / 2.; // sin^2(theta/2) = (1 - cos(theta)) / 2
310  double c2_half = 1. - s2_half; // cos^2(theta / 2) = 1 - sin^2(theta / 2)
311 
312  // Calculate a nonzero cross section only for incident (anti)electrons
313  // or (anti)neutrinos
314  int abs_probe_pdg = std::abs(probe_pdg);
315 
316  /// \todo Add any needed changes for positron projectiles
317  if ( probe_pdg == genie::kPdgElectron ) {
318  // (e,e') differential cross section
319 
320  double q2 = std::pow(q0, 2) - q_mag2;
321 
322  double mott = std::pow(
323  genie::constants::kAem / (2. * E_probe * s2_half), 2) * c2_half;
324 
325  // Longitudinal part
326  double l_part = std::pow(q2 / q_mag2, 2) * entry.W00;
327 
328  // Transverse part
329  double t_part = ( (2. * s2_half / c2_half) - (q2 / q_mag2) ) * entry.Wxx;
330 
331  xsec = (2. * genie::constants::kPi) * mott * (l_part + t_part);
332  }
333  else if ( abs_probe_pdg == genie::kPdgNuE
334  || abs_probe_pdg == genie::kPdgNuMu
335  || abs_probe_pdg == genie::kPdgNuTau )
336  {
337  // Simplify the expressions below by pre-computing the structure functions
338  double w1 = this->W1(q0, q_mag, entry);
339  double w2 = this->W2(q0, q_mag, entry);
340  double w4 = this->W4(q0, q_mag, entry);
341  double w5 = this->W5(q0, q_mag, entry);
342 
343  // Flip the sign of the terms proportional to W3 for antineutrinos
344  double w3 = this->W3(q0, q_mag, entry);
345  if (probe_pdg < 0) w3 *= -1;
346 
347  double part_1 = (2. * w1 * s2_half) + (w2 * c2_half)
348  - (w3 * (E_probe + El) * s2_half);
349 
350  double part_2 = (w1 * cos_l) - (w2/2. * cos_l)
351  + (w3/2. * (El + k_final - (E_probe + El)*cos_l))
352  + (w4/2. * (std::pow(ml, 2)*cos_l + 2*El*(El + k_final)*s2_half))
353  - (w5/2. * (El + k_final));
354 
355  double all_terms = part_1 + std::pow(ml, 2)
356  * part_2 / (El * (El + k_final));
357 
358  xsec = (2. / genie::constants::kPi) * k_final * El
359  * genie::constants::kGF2 * all_terms;
360  }
361 
362  return xsec;
363 }
364 
366  const Interaction* interaction, double Q_value) const
367 {
368  // Don't do anything if you've been handed a nullptr
369  if ( !interaction ) return 0.;
370 
371  int probe_pdg = interaction->InitState().ProbePdg();
372  double E_probe = interaction->InitState().ProbeE(kRfLab);
373  double m_probe = interaction->InitState().Probe()->Mass();
374  double Tl = interaction->Kine().GetKV(kKVTl);
375  double cos_l = interaction->Kine().GetKV(kKVctl);
376  double ml = interaction->FSPrimLepton()->Mass();
377 
378  return dSigma_dT_dCosTheta_rosenbluth(probe_pdg, E_probe, m_probe, Tl, cos_l, ml,
379  Q_value);
380 }
381 
383  double E_probe, double m_probe, double Tl, double cos_l, double ml,
384  double Q_value) const
385 {
386  // dSigma_dT_dCosTheta in GeV^(-3)
387  double xsec = 0.;
388 
389  /// \todo Add check if in grid. If not, return 0
390 
391  // Final state lepton total energy
392  double El = Tl + ml;
393 
394  // Energy transfer (uncorrected)
395  double q0 = E_probe - El;
396 
397  // The corrected energy transfer takes into account the binding
398  // energy of the struck nucleon(s)
399  double q0_corrected = q0 - Q_value;
400 
401  // Magnitude of the initial state lepton 3-momentum
402  double k_initial = real_sqrt( std::pow(E_probe, 2) - std::pow(m_probe, 2) );
403 
404  // Magnitude of the final state lepton 3-momentum
405  double k_final = real_sqrt( std::pow(Tl, 2) + 2*ml*Tl );
406 
407  // Magnitude of the 3-momentum transfer
408  double q_mag2 = std::pow(k_initial, 2) + std::pow(k_final, 2)
409  - 2.*k_initial*k_final*cos_l;
410  double q_mag = real_sqrt( q_mag2 );
411 
412  // Find the appropriate values of the hadron tensor elements for the
413  // given combination of q0_corrected and q_mag
414  TableEntry entry = fGrid.interpolate(q0_corrected, q_mag);
415 
416  // The half-angle formulas come in handy here. See, e.g.,
417  // http://mathworld.wolfram.com/Half-AngleFormulas.html
418  double s2_half = (1. - cos_l) / 2.; // sin^2(theta/2) = (1 - cos(theta)) / 2
419  double c2_half = 1. - s2_half; // cos^2(theta / 2) = 1 - sin^2(theta / 2)
420 
421  // Calculate a nonzero cross section only for incident (anti)electrons
422  // or (anti)neutrinos
423  int abs_probe_pdg = std::abs(probe_pdg);
424 
425  // This Q^2 is defined as negative (Q^2<0)
426  double q2 = std::pow(q0, 2) - q_mag2;
427 
428  /// \todo Add any needed changes for positron projectiles
429  if ( probe_pdg == genie::kPdgElectron ) {
430  // (e,e') differential cross section
431 
432  // double mott = std::pow(
433  // genie::constants::kAem / (2. * E_probe * s2_half), 2) * c2_half;
434  // the previous expression was an approximation in the case of ml-->0 (UR limit).
435  // To be more accurate, SIGMA_MOTT SHOULD BE EXPRESSED AS:
436  double mott = std::pow(2.*El*genie::constants::kAem / q2, 2);
437  // ALTHOUGH THE DIFFERENCES SHOULD BE VERY SMALL OR NEGLIGIBLE
438 
439  // Longitudinal part
440  double l_part = std::pow(q2 / q_mag2, 2) * entry.W00 * c2_half;
441 
442  // Transverse part
443  double t_part = ( 2.*s2_half - (q2*c2_half / q_mag2) ) * entry.Wxx;
444 
445  // This is the double diff cross section dsigma/dOmega/domega==dsigma/dOmega/dEl
446 
447  //xsec = (2. * genie::constants::kPi) * mott * (l_part + t_part);
448 
449  //Bug in HT means we Wxx is off by factor of 2
450  xsec = (2. * genie::constants::kPi) * mott * (l_part + t_part/2.);
451  }
452  else if ( abs_probe_pdg == genie::kPdgNuE
453  || abs_probe_pdg == genie::kPdgNuMu
454  || abs_probe_pdg == genie::kPdgNuTau )
455  {
456 
457  // sigma_0 (sig0) for neutrinos. Similar to sigma_mott for electrons
458  double v0=4.*El*E_probe+q2; // global factor v_0==cos^2(\tilde\theta/2)
459 
460  // CKM matrix element connecting the up and down quarks
462  ->CommonList("Param", "CKM");
463  double Vud = temp_reg->GetDouble( "CKM-Vud" );
464 
465  double Vud2 = std::pow(Vud, 2);
466  double sig0 = genie::constants::kGF2 * Vud2 * v0 * k_final
467  / (4. * genie::constants::kPi * E_probe);
468 
469  // Definition of dimensionless factors
470  double xdelta=ml/sqrt(-q2); // delta
471  double xrho=-q2/q_mag2;
472  double xrhop=q_mag/(El+E_probe);
473  double tan2th2=-q2/v0;
474 
475  // Definition of the lepton kinematic factors, V_K, related to the lepton tensor
476  double VCC=1-std::pow(xdelta, 2)*tan2th2;
477  double VCL=q0/q_mag+tan2th2*std::pow(xdelta, 2)/xrhop;
478  double VLL=std::pow(q0, 2)/q_mag2+std::pow(xdelta, 2)*tan2th2*(1.+2.*q0/q_mag/xrhop+xrho*std::pow(xdelta, 2));
479  double VT=xrho/2.+tan2th2-tan2th2*std::pow(xdelta, 2)/xrhop*(q0/q_mag+std::pow(xdelta, 2)*xrho*xrhop/2.);
480  double VTP=tan2th2/xrhop*(1-q0/q_mag*xrhop*std::pow(xdelta, 2));
481 
482  // Definition of the hadron (nuclear) response functions, R_K, related to the hadron (nuclear) tensor, Wij.
483 
484  double RCC=entry.W00;
485  double RCL=-entry.ReW0z; // RCL must be always negative
486  double RLL=entry.Wzz;
487  double RT=2.*entry.Wxx;
488  double RTP=-entry.ImWxy; // RTP must be positive for nu and negative for nubar
489  if (probe_pdg < 0) RTP *= -1; // THIS IS FOR NUBAR
490 
491 
492  // Determination of the double differential cross section: dsigma/dcostheta_ldEl.
493  // In order to calculate dsigma/dcostheta_ldp_l, the c.s. must be multiplied by
494  // k_final/El
495  xsec= sig0*(VCC*RCC+2.*VCL*RCL+VLL*RLL+VT*RT+2.*VTP*RTP);
496 
497 
498  // This should never happen using the full SuSAv2-MEC hadron tensors
499  // but can trigger when using the tensors from the parameterisation
500  if ( xsec < 0 ) {
501  xsec=0.;
502  }
503  }
504 
505  return xsec;
506 }
const int kPdgNuE
Definition: PDGCodes.h:28
virtual std::complex< double > tz(double q0, double q_mag) const
The tensor element .
#define pERROR
Definition: Messenger.h:59
const int kPdgNuMu
Definition: PDGCodes.h:30
virtual std::complex< double > tt(double q0, double q_mag) const
The tensor element .
TParticlePDG * Probe(void) const
RgDbl GetDouble(RgKey key) const
Definition: Registry.cxx:474
const int kPdgElectron
Definition: PDGCodes.h:35
Registry * CommonList(const string &file_id, const string &set_name) const
Summary information for an interaction.
Definition: Interaction.h:56
virtual std::complex< double > zz(double q0, double q_mag) const
The tensor element .
virtual double dSigma_dT_dCosTheta(const Interaction *interaction, double Q_value) const
virtual double W2(double q0, double q_mag, double Mi) const
virtual double W4(double q0, double q_mag, double Mi) const
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
virtual double W1(double q0, double q_mag, double Mi) const
The structure function .
int A() const
Mass number of the target nucleus.
const Kinematics & Kine(void) const
Definition: Interaction.h:71
int ProbePdg(void) const
Definition: InitialState.h:64
double GetKV(KineVar_t kv) const
Definition: Kinematics.cxx:323
void read1DGridValues(int num_points, int flag, std::ifstream &in_file, std::vector< double > &vec_to_fill)
virtual double W6(double q0, double q_mag, double Mi) const
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
virtual std::complex< double > xx(double q0, double q_mag) const
The tensor element .
void set_pdg(int pdg)
Set the target nucleus PDG code.
virtual double W5(double q0, double q_mag, double Mi) const
virtual double dSigma_dT_dCosTheta_rosenbluth(const Interaction *interaction, double Q_value) const
const int kPdgNuTau
Definition: PDGCodes.h:32
virtual std::complex< double > xy(double q0, double q_mag) const
The tensor element .
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:65
int IonPdgCode(int A, int Z)
Definition: PDGUtils.cxx:71
int Z() const
Atomic number of the target nucleus.
const InitialState & InitState(void) const
Definition: Interaction.h:69
double ProbeE(RefFrame_t rf) const
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
TabulatedLabFrameHadronTensor(const std::string &table_file_name)
static AlgConfigPool * Instance()
virtual double W3(double q0, double q_mag, double Mi) const